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Preface 

This  Is  the  final  report  on  activities  and  research  performed 
under  Contract  DAAD07-78-C-0063.  Most  of  the  work  done  under  tfca$,  con- 
tract has  been  directed  toward  searches  for  ways  in  which  existing  ASir>— 
aerosol  IR  transmittance  modeling  codes  could  be  made  computationally 
more  efficient  in  order  to  enhance  their  utility  in  applications 
involving  time  dependent  variations  in  aerosol  properties.  In  the  course 
of  those  studies,  NMSU  personnel  have  examined  both  rigorous  and 
approximate  methods  for  calculating  both  single  and  multiple  scattering 
extinction  and  scattering  effects,  and  have  devised  a new  version  of  the 
NMSU/ASL  code  AGAUS  which  permits  the  user  to  "trade-off"  overall 
accuracy  for  significantly  faster  program  execution.  The  latter  code 
has  also  been  modified  in  such  a way  that  it  can  be  easily  used  for 
generating  data  which  can  be  used  directly  by  the  ASL  Smoke  Obscuration 
Model. 


A serious  effort  has  also  been  made  to  generate  Improved  documen- 
tation for  program  AGAUS.  The  new  documentation  will  probably  also 
enhance  the  general  utility  of  the  code.. 

Other  minor  efforts  under  the  contract  have  resulted  in  modifica- 
tion of  the  NMSU/ASL  code  ATRANX  [1]  to  permit  users  to  choose  an 
additional  line-profile  (the  Van  Vleck-Weisskopf  line  shape) , to 
reorganize  certain  portions  of  the  code  to  decrease  computation  times, 
to  permit  users  of  the  aerosol  sub-code  to  allow  for  changes  in  particle 
sizes  with  changes  in  relative  humidity,  and  the  integration  of  a more 
sophisticated  Galatry  code  for  the  generalized  Voigt  profile  case. 

The  work  summarized  above  will  be  presented  in  more  detail  in  the 
main  body  of  this  report. 
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Chapter  1 


SINGLE-SCATTERING  INVESTIGATIONS 
Introduction 


This  chapter  presents  a summary  of  work  characterized  as  single 
scattering  studies  which  was  carried  out  under  the  contract.  These 
studies  have  led  to  the  development  of  two  new  versions  of  the  ASL/NMSU 
single-scattering  code  (AGAUS),  comparisons  of  several  different  Mie 
codes  vis  a vis  computational  speed  and  accuracy,  the  development  of  a 
new  analytic  approximation  for  computation  of  angular  intensity 
distributions  for  polydisperse  aerosols,  and  some  extensions  of  earlier 
studies  of  the  way  in  which  relative  humidity  may  affect  infrared 
transmittance  predictions.  In  addition,  both  newly  developed  versions  of 
the  ASL/NMSU  single-scattering  code  AGAUS  have  been  subjected  to  further 
evaluations  with  respect  to  the  results  generated  by  other  codes,  and 
extensive  documentation  has  been  prepared  as  an  appendix  to  this  report. 


1.1  PROGRAM  AGAUS 9 


1.1.1  Introduction:  Program  AGAUS9  is  an  extended  version  of  earlier 

programs  developed  by  NMSU  and  ASL  code  PCAUSS.  AGAUS 9 has  the 
capability  of  generating  "scattering  fractions"  with  a normalization 
consistent  with  that  used  in  the  ASL  Smoke  Obscuration  Model  (ASLSOM) , 
while  retaining  the  option  of  generating  phase-function  expansion 
coefficients  which  may  be  desired  for  subsequent  usage  in  the  ASL  multiple 
scattering  code  STAR04  (and  its  derivative  codes)  or  the  NMSU/ASL  aerosol 
thermal  emission  code  CLEM70.  Additional  features  to  be  found  in  AGAUS9 
(as  compared  to  PGAUSS  or  earlier  AGAUS  versions)  are: 

(1)  Provision  for  automatic  cycling  over  a range  of  wavelengths 
and  averaging  of  results  over  such  a range  of  wavelengths. 

(2)  Provision  for  treating  changes  in  hygroscopic  particle  sizes 
which  may  occur  with  variations  in  relative  humidity. 

(3)  Inclusion  of  an  internal  subroutine  to  provide  automatic 
look-up  and/or  interpolation  of  optical  constants  for  liquid  water  at 
wavelengths  between  0.35  ym  and  200  ym. 

(4)  The  addition  of  two  new  size-distribution  models: 

A)  User-supplied  Blmodal  parameters, 

B)  Marshall-Palmer  "rain"  model. 

(5)  Conversion  of  most  input  data  to  common  data  formats. 

(6)  Replacement  of  the  undocumented  ASL  "Fog  model"  by  a special 
version  of  the  "modified  gamma"  distribution  whose  control  parameters 
include  "liquid  water  content"  rather  than  droplet  number  densities. 
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(7)  Options  whereby  users  can  either  define  particle  number 
densities  or  let  them  be  calculated  from  mass  density  and  mass  concen- 
tration input  parameters. 

(8)  Addition  of  internal  checks  to  warn  users  if  certain  computed 
quantities  seem  to  be  failing  to  converge,  and  extended  error  messages 
associated  with  failures  in  the  Mle  routine. 


1.1.2  Relationships  between  Scattering  Fractions  and  Phase  Functions 

Among  the  major  objectives  of  the  work  done  under  this  contract 
was  the  conversion  of  the  ASL/NMSU  single  scattering  code  AGAUS  into  a 
form  which  produced  the  types  of  data  required  by  the  ASL  Smoke  Obscuration 
Model  (ASLSOM)  while  retaining  the  variety  of  options  previously 
developed  for  AGAUS.  That  conversion  required  that  AGAUS  be  given  the 
additional  capabilities  of  predicting  (a)  attenuation  coefficients  or 
cross-sections  in  units  of  square  meters  per  milligram  of  aerosol 
material,  and  (b)  so-called  "scattering  fractions"  per  unit  aerosol  mass. 
Unfortunately,  the  available  documentation  for  the  Mle  code  used  in 
ASLSOM' s parent  code  was  not  very  clear  as  to  the  precise  relationships 
between  "scattering  fractions"  and  the  customary  quantities  used  in  Mie 
theory.  That  fact  made  the  conversion  of  AGAUS  a great  deal  more  tedious 
and  time  consuming  than  it  might  otherwise  have  been.  In  order  to 
attempt  to  avoid  future  possible  confusions  on  those  relationships,  they 
will  be  summarized  below  before  proceeding  to  a discussion  of  the 
evaluation  of  the  versions  of  AGAUS  developed  under  this  contract. 

The  so-called  "phase  functions",  denoted  here  by  Pj(vO  and 

"scattering  fractions"  S ^ ( 0)  are  both  derived  from  the  average  scattered 

intensities  (I^+l2)/2  defined  in  Mie  theory.  The  two  quantities  are 

related  to  one  another  for  calculations  associated  with  a single  wave- 
length by  a simple  multiplicative  factor,  but  they  have  different 
interpretations  and  applications. 

Let  I(a,m,k,9)  be  the  average  of  the  intensities  1^(6)  and  i^(9)  for 

scattering  at  angle  9 from  a sphere  whose  Mle-slze  parameter  is 
a • 2irr/l,  and  whose  complex  index  of  refraction  is  n « m-ik.  Furthermore, 
let  n(r)dr  be  the  relative  number  of  aerosol  particles  with  radii  between 
r and  (r+dr),  and  let  Q (a,m,k)  and  Q (a,m,k)  respectively  be  the 

6XC  SCa 

total  extinction  and  scattering  efficiency  factors  as  defined  by  van  de 
Hulst  [2]. 

Now,  define  the  quantities  I,  cext  and  Csca  as  follows: 
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1(9) 


1 f 


N„ 


I(a,m,k,9)n(r)dr, 
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(1.1-1) 


(1.1-2) 
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where 


Nt  - j n(r)dr. 
r-o 


(1.1-4) 


The  "phase-function"  as  defined  in  program(s)  AGAUS  is  given  by 


Pf(9) 


1(9)  , 
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ext 


(1.1-5) 


and  has  the  property  that 
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P,(9)d<t>  d(cos9)  =■  4it-^  = 4tt5  , (1.1-6) 
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where  5 is  called  "the  albedo  for  single  scattering." 
o 

The  "scattering  fractions"  S^(9)  are  defined  by 


Sf (9)  - O^)*  Nt  1(9), 


(1.1-7) 


and  S^(9)  has  the  property 
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2ir  n 

j [ Sf<6)  d*  d(cos8)  - C9ca-NT  • (C9Ca)tottl,  (1.1-8) 

$■0  0*o 

wherein  (^gca^totaj  the  total  scattering  cross-section  per  unit 
mass  of  aerosol  material. 

Thus,  it  will  be  seen  that 
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The  ASL  Smoke  Obscuration  Model  (ASLSOM)  is  coded  with  the  assump- 
tion that  (Cgca)totai  as  derived  from  (9)  will  have  units  of  square 

meters  per  milligram.  ASL/NMSU  programs  AGAUS9  and  AGAUS10,  (see  below) 
on  the  other  hand  are  coded  in  what  are  basically  CGS  units.  Conversion 
from  1(9)  to  S^(9)  in  the  ASLSOM  normalization  therefore  requires  a 
unit-conversion  factor  for  X2  from  la  to  m,  and  for  N_  from  cm-3  to  m-3, 
i.e.,  l 

[X(ym)]2  N_ ( cm- 3 ) x (10~*  S-)*  x 10s  —3  . (1.1-10) 

1 j 1TT1  |Q 


Thus,  a factor  10  6 is  required  to  convert  the  scattering  fractions 
from  the  internal  units  of  NMSU/ASL  program  AGAUS  to  the  units  expected 
by  ASLSOM. 


1.1.3  Evaluation  of  AGAUS9:  Once  the  relationships  between  the  Mie 
intensities,  i^  and  i 2 described  above  were  clearly  elucidated,  it 

seemed  that  it  would  be  a straightforward  task  to  convert  AGAUS  from  the 
computation  of  phase  functions  to  the  calculation  of  scattering  fractions. 
When  the  appropriate  conversions  were  completed,  the  new  version,  AGAUS9, 
was  run  using,  as  nearly  as  could  be  determined,  the  same  aerosol  model 
for  which  results  were  found  in  the  SOM  documentation.*  The  particular 


The  label  SOM  is  used  to  refer  to  the  document  entitled  "The  Effectiveness 
of  Obscuring  Smokes",  by  Johnson,  Forney  and  Dolce,  an  unpublished  descrip- 
tion of  the  JTCG/ME  smoke  obscuration  model  which  was  supplied  to  NMSU 
by  R. B.  Gomez  of  ASL. 


model  used  was  a log  normal  distribution  for  WP  smoke  (r  - 0.37  um, 

0 ■ 1.54  Urn),  200  particle  radii  between  minimum  and  maximum  values 
of  0.005  and  1.0  um,  respectively,  a mass  density  of  1.87  gm/cc,  a 
particle  number  density  of  1.276345  x 103  cm-3,  and  complex  Index 
of  refraction  m ■ 1.43  - 01.  Computations  were  performed  at  7 
wavelengths  (0.40,  0.45,  0.50,  0.55,  0.60,  0.65  and  0.70  um) , and  the 
scattering  fractions  were  arithmetically  averaged  over  wavelength. 

Some  of  the  results  are  presented  in  Table  1.1.1.  Only  a cursory 
examination  of  that  table  is  needed  to  observe  that  the  AGAUS9  results 
were  not  in  very  good  agreement  with  those  found  in  the  SOM  documenta- 
tion. This  comparison  was  something  of  a shock,  because  the  basic 
AGAUS  code  had  previously  been  found  to  give  much  better  agreement  with 
other  Mie  codes  than  is  seen  In  Table  1.1.1. 

Because  thorough  reviews  of  the  program  listings  failed  to  shed 
any  light  on  the  reasons  for  such  large  apparent  discrepancies  it  was 
decided  that  comparisons  of  computations  made  at  a single  wavelength 
(rather  than  after  averaging  over  several  wavelengths)  were  highly 
desirable.  The  problem  there  was  that  the  SOM  documentation  did  not 
contain  any  single  wavelength  results.  Since  that  avenue  for  testing 
of  AGAUS 9 was  closed,  some  other  approach  was  needed,  and  was  found 
through  the  courtesy  of  ASL  which  supplied  a copy  of  the  Radiation 
Research  Associates'  code  MIE-2. 

The  aerosol  model  used  to  generate  Table  1.1.1  was  then  passed 
through  both  AGAUS9  and  the  ASL  supplied  MIE-2  code  using  a wavelength 
of  0.40  ym  only.  The  results  of  those  two  runs  are  shown  in  Table  1.1.2.* 

Examination  of  Table  1.1.2  reveals  that  the  AGAUS9  results  and 
the  MIE-2  results  are  in  much  better  agreement  than  those  seen  in 
Table  1.1.1.  Whereas  Table  1.1.1  shows  discrepancies  even  in  the  first 
digit  at  many  scattering  angles,  Table  1.1.2  shows  that  AGAUS9  and 
MIE-2  agree  to  at  least  five  and  often  six  significant  digits — the 
nominal  accuracy  which  is  used  in  terminating  Mie  series  calculations 
in  AGAUS 9.  It  appears,  therefore,  that  the  disagreements  found  in 
Table  1.1.1  should  be  interpreted  as  casting  doubts  on  the  results  quoted 
in  the  SOM  documentation  rather  on  AGAUS9.** 


* 

The  results  of  the  ASL  supplied  version  of  MIE-2  have  been  multiplied 
by  a factor  of  102  to  offset  different  normalizations  used  in  SOM  and 
AGAUS9  and  ASL's  version  of  MIE-2. 

kk 

It  is  possible,  but  not  demonstrated  at  this  point,  that  the  discrepancies 
might  be  traced  to  the  use  of  a "single-precision"  version  of  MIE-2  in 
generating  the  data  listed  in  the  SOM  documentation.  Both  AGAUS9  and  the 
ASL  supplied  version  of  MIE-2  used  in  generating  table  1.1.2  are  "double- 
precision" codes. 


Table  1.1.1  Comparison  of  Scattering  Fractions  from 
AGAUS9  and  SOM  Report  as  Averaged  for  a 
Seven  Wavelength  Run. 


Scattering 

Angle 

Scattering  Fractions 
(NMSU) 

Scattering  Fractions 
SOM  report 

0 

0.5840074-02 

0.530822-02 

10* 

0.3385669-02 

0.294243-02 

20* 

0.1236144-02 

0.961893-03 

30* 

0.5525204-03 

0.419528-03 

40* 

0.2882904-03 

0.254840-03 

50* 

0.1637546-03 

0.176892-03 

60* 

0.9840574-04 

0.114931-03 

70* 

0.6277324-04 

0.706178-04 

80* 

0.4262488-04 

0.445657-04 

90* 

0.3097364-04 

0.309732-04 

100* 

0.2429897-04 

0.249208-04 

110* 

0.2086846-04 

0.237614-04 

120* 

0.2023009-04 

0.265251-04 

130* 

0.2267312-04 

0.330425-04 

140* 

0.3014750-04 

0.446135-04 

150° 

0.4592933-04 

0.637035-04 

160* 

0.6781825-04 

0.883152-04 

170* 

0.4512012-04 

0.678717-04 

180* 

0.6003613-04 

0.836141-04 

T 
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Table  1.1.2 

Comparison  of  Results  from  Runs  of 
ACAUS9  and  MIE2  at  a Single  Wavelength 

Scattering 

Angle 

Scattering  Fractions 
(NMSU) 

I*  (RRA-MIE2 
avg 

0* 

8.215439-03 

8.215442-03 

10 

3.341243-03 

3.341238-03 

20 

8.565565-04 

8.565530-04 

30 

4.126651-04 

4. 126640-04 

40 

2.452253-04 

2.452247-04 

50 

1.503114-04 

1.503111-04 

60 

9.399925-05 

9.399906-05 

70 

6.098411-05 

6.098398-05 

80 

4.167272-05 

4.167267-05 

90 

3.019892-05 

3.019886-05 

100 

2.333270-05 

2.333268-05 

110 

1.957286-05 

1.957282-05 

120 

1.854355-05 

1.854351-05 

130 

2.096545-05 

2.096542-05 

140 

2.942918-05 

2.942910-05 

150 

4.947500-05 

4.947498-05 

160 

8.643879-05 

8.643865-05 

170 

6.208783-05 

6.208780-05 

180 

7.644555-05 

7.644557-05 

2 

The  RRA-MIE-2  results  have  been  multiplied  by  10  to  normalise  them 
In  the  same  way  as  found  in  the  SOM  document. 
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To  avoid  drawing  erroneous  conclusions  on  the  cross-agreements 
between  results  produced  by  AGAUS9  and  MIE-2  on  the  basis  of  a single 
aerosol  model,  additional  comparisons  were  made  using  a quite  different 
aerosol  model.  The  model  chosen  for  this  "test"  was  a "modified"  gamma 
distribution: 


f(r)  - ArQ  exp(-Br^),  [for  MIE-2)  (1.1-11) 

with  A ■ 1.0396x10®,  a ■ 7.5,  B • 28.333333  and  y “ 1.0.  The  input 
data  used  by  AGAUS9  are  the  particle  number  density  NT  and  mode  radius 
rc,  rather  than  A and  B,  but  the  various  quantities  are  related  to  one 
another  through: 

1/Y 

CL  ' ' * v ' Ct¥\ 

rc  - <^J)  and  NT  - AB  y ?(—)-  (1.1-12) 

For  this  comparison  the  minimum  and  maximum  particle  radii  were  taken 
to  be  0.01  pm  and  1.5  urn,  respectively,  and  200  Individual  values  of 
particle  radii  were  used.  The  run  was  made  at  a wavelength  of  0.6  um 
using  m ■ 1.53  - 0.0061.  Some  results  of  the  first  runs  of  this  tvpe 
will  be  found  in  Table  1.1.3. 

Contrary  to  the  excellent  agreement  found  between  MIE-2  and  AGAUS9, 
the  data  found  in  Table  1.1.3  did  not  agree  as  well  as  had  been  expected 
being  good  only  to  about  4 digits.  Further  analysis  indicated  that 
the  most  probable  source  of  these  disagreements  lay  in  the  fact  that 
AGAUS9  and  MIE-2  did  not  choose  the  values  of  particle  radii  actually 
used  in  the  calculation  in  the  same  way.  Slight  alterations  to  remove 
these  differences  were  then  made  in  AGAUS9,  and  the  calculations  were 
repeated,  yielding  the  data  presented  in  Table  1.1.4. 

In  Table  1.1.4  discrepancies  between  the  AGAUS9  and  MIE-2  results 
again  appear  only  in  the  fifth  or  sixth  significant  digit. 

The  above  results  Indicate  once  again  that  code  ACAUS9  is 
quite  capable  of  yielding  results  which  are  as  reliable  as  those 
of  MIE-2.  The  changes  in  AGAUS9  needed  to  bring  about  the  level  of 
agreement  seen  in  Table  1.1.4  involved  only  the  choice  of  the  values  of 
the  200  radii  used,  and  illustrates  that  the  method  of  choosing  the  radii 
can  have  a significant  effect  on  the  values  of  the  scattering  fractions 
even  when  a relatively  large  number  of  radii  are  used.  That  matter  will 
be  discussed  further  in  a subsequent  section  of  this  report. 

At  this  point,  AGAUS9  has  been  shown  to  be  as  valid  as  MIE-2  in 
the  compx’.tation  of  scattering  fractions,  so  that  the  conversion  of 
AGAUS  into  an  ASLSOM  compatible  form  has  been  successfully  accomplished. 


Table  1.1.3  Comparison  of  Scattering  Fractions 

Calculated  by  MIE-2  and  AGAUS9  for  a 
Modified  Gamma  Distribution 


Scattering 

Angle 

Scattering  Fractions 
CNMSU) 

I*  (RRA-MIE2) 
avg' 

0° 

9.254421-08 

9.254768-08 

10 

7.620680-08 

7.620913-08 

20 

4.470319-08 

4.470374-08 

30 

2.156754-08 

2.156747-08 

40 

1.038258-08 

1.038255-08 

50 

5.626184-09 

5.626226-09 

60 

3.415889-09 

3.415935-09 

70 

2.236736-09 

2.236776-09 

80 

1.555144-09 

1.555170-09 

90 

1.149602-09 

1.149602-09 

100 

9.113908-10 

9.114065-10 

110 

7.876338-10 

7.876470-10 

120 

7.656564-10 

7.656688-10 

130 

8.657526-10 

8.657658-10 

140 

1.074210-09 

1.074245-09 

150 

1.193137-09 

1.193200-09 

160 

1.134702-09 

1.134755-09 

170 

1.502157-09 

1.502236-09 

180 

1.969069-09 

1.969203-09 

x 2 

The  RRA-MIE-2  results  have  been  multiplied  by  10  to  normalize  them 
in  the  same  way  as  found  in  the  SOM  document. 
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Table  1.1.4  Comparison  of  Scattering  Fractions  Produced 
by  MIE-2  and  a Special  Version^  of  AGAUS9 
for  a "Modified  Gamma"  Distribution  Model. 


Scattering 

Angle 

Scattering  Fractions 
(NMSU) 

I*  (RRA-MIE2) 
avg 

0* 

9.254772-08 

9.254768-08 

10 

7.620918-08 

7.620913-08 

20 

4.470380-08 

4.470374-08 

30 

2.156751-08 

2.156747-08 

40 

1.038257-08 

1.038255-08 

50 

5.626233-09 

5.626226-09 

60 

3.415938-09 

3.415935-09 

70 

2.336777-09 

2.336776-09 

80 

1.555171-09 

1.555170-09 

90 

1.149623-09 

1.149622-09 

100 

9.114071-10 

9.114065-10 

110 

7.876480-10 

7.876470-10 

120 

7.656697-10 

7.656688-10 

130 

8.657668-10 

8.657658-10 

140 

1.074245-09 

1.074245-09 

150 

1.193200-09 

1.193200-09 

160 

1.134755-09 

1.134755-09 

170 

1.502236-09 

1.502236-09 

180 

1.969203-09 

1.969203-09 

^particle  radii  selected  the  same  way  as  in  MIE-2. 

* 

The  RRA-MIE-2  results  have  been  multiplied  by  10*  to  normalize  them 
in  the  same  way  as  found  in  the  SOM  document. 


1.2  PROGRAM  AGAUS10 


1.2.1  Introduction:  Experience  acquired  over  a period  of  several 

years  with  various  versions  of  the  single  scattering  codes  PGAUSS  and 
AGAUS  has  revealed  that  users  of  such  codes  often  tend  to  be  over- 
conservative in  their  choice  of  the  number  of  radii  at  which  Mie  calcu- 
lations are  performed.  Many  single  scattering  codes  (including  PGAUSS, 
AGAUS9,  and  R.R.A. 's  MIE-2)  require  that  the  user  specify  some  input 
parameter  such  as  the  total  number  of  radii  to  be  used.  Unless  a user 
has  had  a great  deal  of  experience  with  choosing  such  parameters  for 
various  types  of  size  distributions,  there  is  a tendency  to  use  many 
more  radii  than  may  be  needed  to  obtain  results  at  acceptable  levels  of 
accuracy.  Since  the  overall  running  time  of  Mie  codes  may  be  greatly 
reduced  by  reducing  the  number  of  particle  sizes  treated,  it  is  desirable 
to  have  a code  in  which  one  specifies  an  acceptable  level  of  accuracy  and 
which  then  uses  only  as  many  particle  radius  values  as  are  needed  to 
satisfy  that  requirement.  The  concept  of  a Mie  code  written  along  that 
line  is  not  new  and,  in  fact,  ASL  already  has  a code  of  that  nature 
but  it  is  one  which  is  limited  solely  to  log-normal  size  distributions. 

In  view  of  the  fact  that  substantial  savings  in  computer  time  might  be 
obtained,  while  maintaining  a great  deal  of  versatility,  NMSU  has 
developed  a special  version  of  AGAUS9  called  AGAUS10  which  is  based  on 
the  sorts  of  compromises  indicated  above. 

Code  AGAUS10  is  in  many  ways  identical  to  AGAUS9,  but  the  method 
for  choosing  values  for  particle  radius  and  the  scheme  for  integration 
over  a size  distribution  is  quite  different.  Rather  than  specifying 
the  number  of  radii  (NRADI  for  AGAUS9),  the  user  of  AGAUS10  must 
specify  a "convergence"  level  DELTA.  The  quantity  DELTA  represents  the 
minimum  fractional  accuracy  that  is  acceptable  for  certain  results  of 
a run.  In  operation,  AGAUS10  then  runs  through  the  Mie  calculations  and 
integrations  over  size  using  only  a few  particle  radii  and  obtains  first 
estimates  of  extinction  coefficients,  etc.  The  size  intervals  are  then 
reduced  to  one-half  their  former  value  and  new  estimates  are  calculated. 

If  the  new  and  old  estimates  agree  to  within  (DELTAxlOO)  percent,  then 
AGAUS 10  ceases  to  treat  additional  particle  sizes  and  proceeds  to  its 
final  calculations.  If  the  new  estimate  does  not  agree  with  the  "old" 
one  to  within  (DELTAxlOO)  percent,  the  size  interval  is  again  cut  in 
half  and  the  calculations  and  comparisons  are  repeated.  The  cycling 
occurs  repetitively  until  either  the  convergence  criterion  is  satisfied 
or  until  some  pre-coded  maximum  number  of  values  (513  for  AGAUS10)  for 
particle  radius  he”®  been  used. 

It  is  believed  that  AGAUS10  will  represent  a substantial  advancement 
in  ASL's  aerosol  modeling  capabilities  because  it  contains  a large  variety 
of  size  distributions,  the  possibility  of  modeling  mixtures  of  aerosols 
with  different  optical  properties,  treats  changes  in  particle  size  with 
variable  relative  humidity,  is  compatible  with  ASLSOM,  and  will  probably 
provide  substantially  shorter  computation  times  than  most  existing 
aerosol  codes. 


. 
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In  the  next  few  sections  of  this  report  will  be  found  a more 
detailed  description  of  the  theoretical  and  numerical  techniques  which 
distinguish  AGAUS10  from  Its  predecessors,  and  comparisons  between  the 
results  and  running  times  AGAUS9 , AGAUS10  and  the  R.R.A.  code  MIE-2. 


1.2.2  Further  Description  of  the  "Halving"  Methods  of  AGAUS10 

Basic  Concept . A general  Idea  of  how  program  AGAUS10  operates 
can  be  obtained  by  considering  a numerical  method  for  determining  the 
area  under  a curve  g(R)  vs  R such  as  that  sketched  In  figure  1.2.1. 
The  objective  Is  to  evaluate  numerically  an  Integral 


G 


R*» 

g(R)  dR 

R-o 


(1.2-1) 


to  some  desired  agree  of  accuracy  using  the  smallest  number  of  values  of 
R for  the  numerical  calculations.  The  procedure  adopted  In  AGAUS10  Is 
as  follows: 


(1)  An  initial  estimate  of  the  value  of  the  integral  is  made 

using  three  values  of  R labeled  by  Roman  numeral  I and  the  "trapezoidal 
rule". 

(2)  A second  estimate  Gj  is  then  made  using  increments  AR  which 
are  half  as  large  as  those  used  in  getting  G^.  In  getting  G.,,  the  two 

additional  R-values  labeled  II  are  utilized. 

(3)  The  values  of  G^  and  G^  are  compared  to  each  other  by  calculating 
a quantity 


(1.2-2) 


and  comparing  it  to  a pre-set  quantity  A. 

If  6 < A,  then  it  is  assumed  that  G0  is  a "sufficiently"  accurate 

representation  of  G,  and  the  computations  are  terminated.  If,  on  the 
other  hand,  5 > A,  one  proceeds. 

(4)  A third  estimate  G^  of  G is  then  made  by  again  cutting  the 

increments  AR  to  half  its  previous  value.  This  results  in  the  addition 
of  computations  at  the  four  new  R-values  labeled  III  in  the  figure  1.2.1. 
A new  value  of  6 is  then  calculated  from 


.1 

i 

I 


* 


■v 
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|g3-g2| 


(1.2-3) 


and  5 la  again  compared  Co  A.  If  6 Is  still  greater  than  A,  the 
spacing  between  R-values  Is  cut  In  half  once  more,  and  the  "estimation- 
comparison"  process  is  repeated  until  either  (a)  6 < A,  or  (b)  some 
maximum  number  of  R-values  has  been  reached. 

In  AGAUS10,  the  quantity  used  in  the  testing  process  is  the  total 
volume  of  the  aerosol  particles,  namely. 


V - J ttR1  f (R)dR,  (1.2-4) 

In  which  f(R)dR  is  the  relative  number  of  aerosol  particles  whose  radii 
lie  between  R and  R+dR. 

A more  formal  description  of  the  mathematical  methods  used  will  be 
found  in  the  appendices. 


1.2.3  Extensions  to  the  Basic  Concept:  It  may  happen  that  the  function 
g(R)  illustrated  in  figure  1.2.1  will  have  a form  which  is  much  "smoother" 
over  some  range  of  R-values  than  over  other  ranges.  Strict  application 
of  the  halving  procedure  to  such  a situation  could  easily  lead  to  the 
use  of  many  more  R-values  in  one  of  those  ranges  than  are  really  needed 
because  of  the  behavior  of  g(R)  in  other  R-ranges.  Such  a situation  seems 
especially  likely  with  asymmetric  sire-distributions.  One  way  of  attempting 
to  avoid  un-needed  computations  in  various  R-ranges  is  to  split  the 
regions  used  in  halving  loops  into  several  distinctive  "intervals" — say, 

R^  to  Rb,  Rjj  to  R^,,  Rc  to  Rq,  ....  and  to  do  the  halving  operations  and 

convergence  tests  separately  for  each  "interval".  AGAUS10  has  been 
coded  to  permit  separate  halving  computations  and  convergence  testing 
for  up  to  eight  "intervals"  of  that  type.  In  fact,  however,  the  present 
coding  uses  two  radius  "intervals"  (ranges)  for  those  distributions 
f(R)  known  to  show  a single  maximum,  four  "intervals"  for  bi-modal 
distribution  models,  and  just  one  interval  for  other  models. 

The  detailed  mathematical  description  of  the  use  of  multiple 
size- Intervals  will  also  be  found  In  the  appendices. 


1.2.4  Usage  of  Program  AGAUS 10:  AGAUS10  has  been  developed  somewhat 
concurrently  with  AGAUS9,  and  a strenuous  effort  has  been  made  to  make 
their  Input  and  output  data  more  or  less  Identical.  It  is  believed  that 
the  objective  has  generally  been  fulfilled,  but  a few  differences  in 
Input  data  exist.  In  particular,  AGAUS9  requires  the  user  to  supply 
a value  (NRADI)  for  the  number  of  radii  to  be  used  In  a run,  whereas 
AGAUS10  requires  that  the  user  specify  a value  for  the  convergence 
level  A (DELTA).  Those  input  data  are  arranged,  however,  so  that  the 
same  set  of  data  cards  can  be  used  as  Input  to  either  program:  AGAUS9 
will  ignore  any  given  value  of  DELTA,  and  AGAUS10  will  ignore  the  user 
supplied  value  of  NRADI — except  for  the  case  IDSTP  - 0,  in  which  f(R) 
is  to  be  supplied  externally  by  the  user. 

Programs  AGAUS9  and  AGAUS10,  for  the  most  part,  use  the  same 
subroutines  (ANGLE,  GAUS,  GUSET,  MIEGX,  WATER),  or  slight  variants 
(AG9PT3  is  the  same  as  AGXPRT  except  in  title;  the  same  is  true  of 
AG9PRT  and  AGXPRT).  Non- interchangeable  routines  are  the  main  programs, 
subroutines  AG9PT1/AGXPT1  and  routines  AG9PT2 /AGXPT2 . 

Sample  compatible  input  data  decks  for  several  distribution  types 
will  be  found  in  appendix  A. 4. 


1.2.5  Evaluation  of  Program  AGAUS10 

1.2. 5.1  Introduction:  In  the  development  of  the  actual  coding  of 
program  AGAUS10,  continuous  testing  and  comparison  of  results  produced 
by  the  new  codes  and  those  yielded  by  other  single  scattering  codes  was 
carried  out.  The  findings  of  a few  such  comparisons  will  be  given  next. 

The  supplementary  codes  used  for  evaluation  of  AGAUS10  were  the 
NMSU/ASL  code  AGADS9,  and  the  Radiation  Research  Associates  code,  MIE-2, 
which  was  supplied  to  NMSU  by  ASL.  In  order  to  avoid  the  problem  of 
deciding  which  of  the  three  codes  was  "correct"  in  the  event  that  no  two 
agreed  well,  all  codes  were  run  using  aerosol  models  for  which  independent 
results  could  be  found  in  Deirmend j ian' s book  Electromagnetic  Scattering 
on  Spherical  Polydispersions  [3],  The  specific  models  used  for  comparison 
runs  were  Deirmend j ian ' s cloud  models  C.l  and  C.3  at  a wavelength  of 
0.7  ym.  Additional  comparison  runs  were  made  at  10.6  urn  using  an  NMSU 
simulated  hygroscopic  smoke  model  [1]  and  several  values  of  relative 
humidity.  The  latter  comparisons  involved  only  program  AGAUS9  and 
AGAUS10.  The  objectives  of  all  comparisons  were  not  only  to  show  that 
AGAUS10  was  working  properly,  but  also  to  determine  whether  or  not  it 
would,  in  fact,  offer  significant  reductions  in  overall  running  time. 


1.2. 5. 2 Water  Cloud  Model  Comparisons:  Comparison  runs  of  codes 
AGAUS9,  AGAUS10  and  MIE-2  for  cloud  models  C.l  and  C.3  were  made  using 
the  same  range  of  radii  quoted  for  Deirmendj ian' s tables  T.36  and  T.30. 


— 
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In  order  Co  assure  chat  any  conclusions  drawn  about  relative  computation 
times  would  be  valid  in  the  context  of  ASLSOM  usage,  scattering 
fractions  were  calculated  at  2*  increments  between  scattering  angles  of 
0®  and  180®  In  these  runs.  AGAUS10  was  also  run  for  model  C.l  using 
several  different  values  of  the  convergence  level  DELTA. 


Table  1.2. 5.1  presents  some  of  the  results  obtained  for  cloud 
model  C.l.  The  scattering  fraction  entries  for  the  "Deirmendj ian" 
column  were  obtained  by  averaging  his  and  P-,,  and  multiplying  by 

the  appropriate  conversion  factor.  The  AGAUS10  run  listed  used 
DELTA  '■  0.01.  Examination  of  the  tabulated  data  shows  that  none  of  the 
three  runs  yields  exactly  the  same  extinction  coefficient  that  was 
printed  by  Deirmendj ian,  but  all  were  within  0.3S  of  his  results. 
Unfortunately,  Deirmendj ian' s table  T. 36  did  not  include  angles  at  2® 
increments,  so  only  forward  and  backward  scattering  fractions  could  be 
inferred  accurately.  It  will  be  seen  that  all  three  codes  gave 
scattering  fractions  which  agreed  well  in  the  0 ■ 0®  direction,  but 
the  results  are  obviously  not  identical.  The  variations  are,  however, 
quite  a bit  larger  at  other  scattering  angles,  although  the  variations 
are  really  not  terribly  large  when  their  magnitudes  compared  to  chat 
of  the  forward  scattering.  For  this  particular  model,  it  will  be  noted 
that  AGAUSlO's  convergence  criterion  led  to  its  use  of  more  radii  than 
were  used  with  the  other  codes.  It  will  also  be  seen  that  MIE-2 
handled  the  40  radii  faster  than  did  AGAUS9.  MIE-2' s speed  advantage 
over  AGAUS9  probably  results  from  the  fact  that  MIE-2  uses  a non-uniform 
spacing  of  particle  radii — a procedure  which  leads  to  a significantly 
smaller  number  of  Mie  calculations  which  must  be  done  at  large  Mie 
size-parameters . 

It  should  be  noted  that  the  version  of  AGAUS10  used  in  generating 
the  results  shown  in  Table  1.2.5. 1 performed  its  "convergence  test"  on 
the  total  aerosol  volume  and  not  on  the  scattering  fractions.  (Results 
obtained  by  checking  on  both  volume  and  the  scattered  intensity  at 
90®  will  be  given  below.) 


Although  AGAUS10  appears  to  offer  no  computation  speed  advantage 
if  Che  user  KNOWS  how  many  radii  are  "enough",  it  does  have  the  useful 
feature  (not  found  in  either  MIE-2  or  AGAUS9)  of  giving  the  user  some 
idea  of  how  many  radii  were  really  required  to  achieve  some  minimum 
level  of  confidence  in  its  results. 

‘ 1 

By  repeating  runs  of  AGAUS10  and  Model  C.l  using  various  values 
of  DELTA,  it  has  been  possible  to  learn  more  about  the  sensitivities 
of  the  extinction  coefficient  and  scattering  fractions  to  the  number 
of  radii  used  in  a calculation.  Table  1.2. 5. 2 contains  the  results 
of  a few  runs  of  that  type.  The  first  four  entries  represent  runs  in 
which  the  quantity  which  was  tested  for  convergence  was  the  total  volume 
of  the  aerosol  particles. t The  fifth  and  sixth  entries  are  the  results 


Volume  was  used  rather  than,  say,  extinction  cross-section  because  its  RJ 
dependence  appeared  to  make  it  a "more  sensitive"  test  quantity. 


X 


I 


Table  1 2.5.1  Comparisons  of  Results  and  CPU  Times 
Single  Scattering  Codes  - Cloud  Mode 


Table  1.2.5. 2 Comparison  of  AGAUS10  Results  Obtained  with  Four 

Convergence  Test  Levels  Using  Cloud  Model  C.l.  (>  * 0.7pm) 
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of  runs  in  which  testing  was  done  on  both  volume  and  the  scattered 
intensity  of  90°;  note  that  the  sixth  entry  represents  a run  at  45° 
increments  rather  than  2°  increments. 

Examination  of  Table  1.2. 5. 2 shows  that  the  value  of  the  extinction 
coefficient  is  not  nearly  as  sensitive  to  the  choice  of  particle  radius 
values  as  are  the  scattering  fractions  at  fairly  large  scattering  angles, 
especially  near  90°.  The  90°  results  show  "oscillations"  of  as  much  as 
15%  as  the  number  of  radii  changes  through  the  set  (26,  50,  66,  98,  130, 
and  513}. 

It  should  also  be  noted  that  the  AGAUS9  and  MIE-2  runs  shown  in 
Table  1.2.5. 1 used  only  40  radii,  while  Deirmendjian  evidently  used  440 
radii.  Even  with  that  vastly  larger  number  of  radii,  Deirmendjian's 
results  are  not  so  different  from  the  present  results  to  warrant,  in  many 
applications,  an  additional  increase  in  computation  time  of  a factor  of 
eleven.  It  was,  in  fact,  observations  of  that  type  gleaned  from  earlier 
usage  of  AGAUS  versions,  that  eventually  led  to  the  decision  to  develop 
AGAUS10. 

A similar  set  of  comparative  results  obtained  for  cloud  model  C.3 
are  presented  In  Table  1.2.5. 3.  For  this  model,  AGAUS9  and  MIE-2  were 
run  with  the  same  number  of  radii  as  was  used  by  Deirmendjian.  This 
example  illustrates  the  way  in  which  AGAUS10  can  offer  definite  decreases 
in  CPU  time  over  codes  requiring  the  user  to  make  an  a priori  decision 
as  to  the  number  of  radii  needed.  Using  only  66  radii,  AGAUS10  produced 
an  extinction  coefficient  to  within  better  than  1%  of  Deirmendjian's 
value  (with  which  AGAUS9  and  MIE-2  results  are  identical  to  the  four 
digits  given  by  Deirmendjian).  AGAUS10,  (with  DELTA  = 0.01),  however, 
ran  nearly  four  times  faster  than  either  AGAUS9  or  MIE-2.  Furthermore, 
the  AGAUS10  run  at  DELTA  = 0.01  yielded  scattering  fractions  which,  at 
worst,  differed  from  those  found  by  AGAUS9  and  MIE-2  by  the  order  of  10%. 

Cloud  models  such  as  C.l  and  C.3  used  at  short  wavelengths  (-  visible) 
are  severe  tests  of  single  scattering  codes  because  they  typically 
involve  some  rather  large  Mie  size-parameters  (a  ~ 110  for  the  C.l  model). 
Calculations  performed  at  infrared  wavelengths  for  smaller  particles 
than  found  in  typical  clouds  are  not  quite  so  demanding.  That  fact  is 
illustrated  partially,  by  Tables  1.2. 5. 4 and  1.2. 5. 5,  which  were  compiled 
using  an  NMSU  hypothetical  model  for  a hygroscopic  smoke  having  the 
particle  size-distribution  of  white  phosphorous  smoke.  The  model  used 
here  was  denoted  as  model  A'  in  an  earlier  document  [ll.  Table  1.2. 5. 4 
shows  extinction  coefficients  obtained  from  AGAUS9  (using  100  radii), 
and  AGAUS10  (using  DELTA  * 0.01)  at  seven  different  values  of  relative 
humidity.  It  also  presents  the  results  of  a run  of  a special  version 
of  AGAUS10  (labeled  AGAUSXL)  in  which  the  RRA  Mie-routine  DOWN42  was 
used  instead  of  routine  MIEGX.  It  will  be  seen  that  AGAUS10  (and 
AGAUSXL)  needed  only  34  radii  to  reach  extinction  coefficients  within 
0.2%  of  those  found  by  AGAUS9,  but  the  AGAUS9  run  used  340%  more  computer 


Table  1.2. 5. 3 Comparisons  of  Results  and  CPU  Times  for  Three 

Single  Scattering  Codes  - Cloud  Model  C.3.(A  = 0.7cm). 
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Table  1.2. 5. 4 Comparisons  of  Results  from  Programs  AGAUS9 
and  AGAUS10  for  NMSU  Smoke  Model  A'  at 
X ■ 1(L6  um.  (Relative  humidity  ■ OX) 


Extinction  Coefficients  (per  km) 


Relative  Humidity  AGAUS9 
(percent) (100  radii) 


AGAUS  10  ACAUSXL 

(34  radii)  (34  radii) 

[Delta  - 0.01] X Difference* 


0 

0.1838 

0.1835 

0.1835 

0.16 

75 

0.2801 

0.2797 

0.2797 

0.18 

80 

0.3020 

0.3016 

0.3016 

0.17 

85 

0.3340 

0.3335 

0.3335 

0.27 

90 

0.3887 

0.3882 

0.3882 

0.15 

95 

0.5350 

0.5343 

0.5343 

0.15 

99 

1.9295 

1.9268 

1.9268 

0.15 

Total  CPU 
Time 

34.34  sec 

9.83  sec 

34.98  sec 

71.4 

* 

NOTE:  per  cent  differences  are: 


| (AGAUS9-AGAUS10)/AGAUS9|  x 100. 


Table  1.2. 5. 5 Comparison  of  AGAUS10  Results  for  Four  Different 
Convergence  Levels  and  a Simulated  Smoke  Model  - 
X - 10.6  um.  (Relative  humidity  ■ Ot) . 


Convergence  Level  No.  of  Radii  Kext 

(DELTA)  used (,km~1) 


.05 

IS 

1.8277 

.01 

34 

1.8353 

.005 

50 

1.8345 

.001 

130 

1.8377 

Convergence  Level 
(DELTA) 

O 

O 

a 

Scattering  Fractions 

0 - 90° 

6 - 180 

.05 

1.5574E-6 

6 . 9366E-7 

1 .2326E-6 

.01 

1 . 5610E-6 

6 . 9555E-7 

1 .2365E-6 

.005 

1 . 5605E-6 

6. 9533E-7 

1 . 2361E-6 

.001 

1.5621E-6 

6.9613E-7 

l . 2376E-6 

23 
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time  than  did  AGAUS10.  ACAUSXL's  results  were  quite  close  to  those  of 
AGAUS10,  but  the  relative  slowness  of  DOWN42  brought  the  total  time  to 
treat  34  radii  up  to  nearly  the  same  time  needed  by  AGAUS9  to  treat  100 
radii. 

Table  1.2.5. 5 demonstrates  how  the  user's  choice  of  DELTA  affected 
runs  of  the  simulated  smoke  model  at  a relative  humidity  of  0%.  For 
that  model,  changes  in  the  number  of  radii  from  18  to  130  (a  factor  of 
seven!)  did  not  change  either  the  extinction  coefficient  or  the  scattering 
fractions  by  more  than  0.3  percent. 


1.2. 5. 3 Discussion:  One  of  the  decisions  which  had  to  be  made  for 
program  AGAUS10  was  that  of  determining  just  which  "quantity"  or  set  of 
quantities  should  be  used  in  the  test  for  convergence  to  within  the 
level  DELTA.  The  one  which  was  finally  chosen  was  the  total  volume  of 
the  aerosol  particles.  The  volume  was  chosen  because,  being  dependent 
on  the  cubes  of  the  radii,  it  might  be  a little  more  conservative  choice 
than  the  extinction  coefficient,  although,  as  seen  in  Tables  1.2. 5.1  and 
1.2. 5. 2,  it  is  a less  conservative  test  quantity  than  the  scattered 
intensity  at,  say,  9 ■ 90°.  Present  experience  is  limited,  but  the 
use  of  the  volume  seems  to  be  a good  compromise  unless  highly  precise 
scattering  fractions  are  needed.  If  the  latter  situation  arises  in  a 
particular  application,  users  of  the  code  should  have  no  difficulty  in 
coding  in  their  own  preferred  test  quantity.  The  quantities  which  are 
available  at  present  are  extinction,  scattering  and  backscattering 
cross-sections  and  the  average  intensities  at  various  angles.  Careful 
study  of  the  source  listing  for  subroutine  AGXPT2  will  reveal  how  such 
changes  can  be  made. 

A major  surprise  which  was  revealed  by  these  studies  was  the  fact 
that  the  code  MIE-2  used  less  computer  time  to  handle  cloud  models  C.l 
and  C. 3 than  AGAUS9  needed  (for  the  same  number  of  radii).  The  reason 
for  the  surprise  was  that  the  basic  Mie  routine  DOWN42  used  by  MIE-2 
was  known  to  be  usually  appreciably  slower  than  routine  MIEGX  (see 
section  1.3).  It  was,  in  fact,  that  difference  in  speed  which  caused 
AGAUSXL  (see  above,  and  Table  1.2. 5. 4)  to  require  more  than  three  times 
as  long  to  treat  the  smoke  model  than  AGAUS10  needed.  The  reason  for 
MIE-2's  evident  computation  time  advantage  over  AGAUS9  has  been  alluded 
to  above  (section  1.1),  and  appears  to  be  the  result  of  the  differences 
in  the  ways  in  which  the  values  of  radii  (at  which  the  Mie  calculations 
are  to  be  done)  are  chosen  in  the  two  codes.  In  MIE-2,  the  interval 
between  successive  values  of  particle  radius  is  doubled  as  each  additional 
value  is  assigned  (at  least  for  the  generalized  Khirgian-Mazin  distri- 
bution), while  AGAUS9  uses  a constant  increment.  The  MIE-2  method 
therefore  uses  far  fewer  large  values  of  radius  than  the  AGAUS9  method 
does.  That  difference,  and  the  fact  that  the  time  required  within  the 
Mie  subroutines  DOWN42  and  MIEGX  increases  rapidly  with  the  size  to 
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wavelength  ratio,  appears  to  account  for  MIE-2's  (small)  running  time 
advantage  over  AGAUS9.  That  difference  In  how  the  radii  are  chosen 
definitely  accounts  for  the  disagreements  between  the  MIE-2  and 
AGAUS9  values  for  the  scattering  fractions  for  cloud  models  C.l  and  C. 3 
(that  disagreement  Is  quite  appreciable  [see  Table  1.2.5. 1]  in  some 
cases).  The  latter  conclusion  has  been  verified  by  direct  runs  with  a 
version  of  AGAUS9  which  uses  the  MIE-2  method  for  choosing  the  radii 
(the  AGAUS9/MIE-2  disagreements  moved  out  to  the  fifth  or  sixth  digit). 
The  running  time  of  the  special  version  of  AGAUS9  was  240  seconds  as 
compared  to  266  seconds  for  MIE-2,  demonstrating  that  the  use  of  fewer 
large  radii  significantly  decreases  running  time.  That  alteration  of 
AGAUS9  has  not  been  made  permanent,  although  the  changes  needed  are  not 
particularly  troublesome,  because  it  is  not  really  clear  that  the  MIE-2 
method  of  doubling  the  interval  does  not  also  carry  along  an  implication 
that  a larger  number  of  radii  (than  needed  with  AGAUS9)  would  be  needed 
to  secure  adequate  sampling  of  larger  particles  sizes. 

The  questions  raised  above  are,  however,  made  somewhat  academic  by 
AGAUS10  for  many  applications  requiring  single  scattering  calculations 
because  of  AGAUSlO's  clear  speed  advantage  in  handling  distributions  for 
which  no  ji  priori  guide  as  to  how  many  radii  will  be  "enough"  is 
available. 


1.2. 5.4  Some  Effects  of  the  Number  of  Size-Intervals  Used  in  AGAUS10 


In  the  early  development  of  AGAUS10,  it  was  tentatively  assumed 
that  computation  times  might  be  substantially  reduced  by  breaking  the 
total  ranges  of  R-values  to  be  treated  into  more  than  one  "size- interval" 
and  then  performing  the  halving  calculations  and  tests  separately 
within  each  size  interval.  That  assumption  was  based  upon  the  belief 
that  the  use  of  a single  interval  could  result  in  unneeded  computations 
in  some  size  ranges  brought  about  by  slow  convergence  in  other  ranges. 
Consequently,  all  types  of  distribution  functions  which  were  known  to 
be  "peaked",  were  split  into  two  size-intervals:  (a)  one  with  radii 
smaller  than  the  one  at  which  the  distribution  was  a maximum,  and  (b)  one 
containing  all  larger  values  of  particle  radii.  In  some  of  the  test 
runs,  however,  it  was  found  that  the  interval  containing  the  smaller 
radii  required  more  halvings  than  the  other  interval,  but  contributed 
only  a small  fraction  of  the  total  extinction  effects — thereby  perhaps 
wasting  computation  time. 

The  effect  of  using  Just  a single  size  interval,  instead  of  two 
intervals,  has  been  briefly  explored  using  a special  version  of  AGAUS10, 
and  Deirmendjian's  cloud  models  C.l  and  C.3  at  A ■ 0.7  pm.  Computations 
were  made  with  several  different  values  for  the  convergence  level  (A). 
Table  1.2. 5. 6 compares  some  of  the  results  found  with  the  special 
version  (AGAUS10S) , and  the  regular  version.  The  scattering  fraction 
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Table  1.2. 5. 6 


Comparison  of  Results  Produced  for  Cloud  Model  C.l 
(A  » 0.7  van)  by  Several  Variants  of  Program  AGAUS10 
and  Different  Convergence  Levels  (A) . 


Scattering  Fractions 


DELTA 

* 

NI 

NRADI 

Kext 

CPU 

“ o* — 

90° 

180fc 

0.05 

a 

2 

26 

17.47 

45.28 

2.4521 

5. 2019E-5 

7.4914E-4 

b 

1 

9 

" 

16.70 

21.06 

2.2353 

6. 9736E-5 

1. 2946E-3 

c 

1 

9 

16.70 

note-d 

2.2353 

6. 9736E-5 

1. 2946E-3 

0.01 

a 

2 

50 

16.72 

82.81 

2.2404 

4. 8011E-5 

8 . 3856E-4 

b 

i 

17 

16.77 

38.72 

2.2438 

7.0311E-5 

1. 0222E-3 

c 

i 

129 

16.74 

note-d 

2.2386 

5. 6969E-5 

9. 1422E-3 

0.005 

a 

2 

66 

16.74 

127.93 

2.2395 

5.6570E-5 

8 . 7430E-4 

b 

1 

17 

16.77 

37.78 

2.2438 

7. 0311E-5 

1. 0222E-3 

0.001 

a 

2 

130 

16.73 

245.13 

2.2370 

5.4473E-5 

8. 2393E-4 

b 

1 

33 

16.76 

75.71 

2.2398 

6.8529E-5 

9. 3237E-4 

"Reference 

Case" 

c 

2 

513 

16.75 

noted-d 

2.2411 

5. 5880E-5 

8.8647E-4 

Deirmendj ian 

1 

440 

16.73 

2.2370 

8.457E-4 

* 

Notes 


a)  "Normal"  AGAUS10;  convergence  tests  on  aerosol  volume;  2 size  interval. 

b)  AGAUS10S;  convergence  tests  on  extinction  coefficient;  1 size  interval. 

c)  AGAUS10S;  convergence  tests  on  extinction  coefficient  and  scattered 
intensity  at  0 ■ 90";  number  of  size  intervals  given  in  column 
labeled  "NI". 

d)  45°  angular  increment  runs;  CPU  times  not  directly  comparable  to 
other  runs  made  at  2®  increments. 
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computations  were  made  at  2*  Increments  in  most  cases  shown  In  the 
table,  and  the  column  labeled  "NI"  explicitly  shows  the  number  of 
size-intervals  used.  Also  shown,  and  labeled  "REF.  CASE",  are  results 
obtained  with  an  AGAUS10  run  which  used  all  513  available  values  of 
particle  radii. 

The  results  seen  in  Table  1.2.5. 6 Indicate  that  the  special 
version  (AGAUS10S)  found  the  convergence  of  the  extinction  coefficient 
to  be  satisfied  with  a substantially  smaller  number  of  particle  sizes 
than  the  "normal"  version  required.  That  fact  also  meant  that 
AGAUS10S  required  only  one-third  to  one-half  as  much  computation  time 
to  satisfy  the  different  convergence  levels  as  AGAUS10  needed.  It 
thus  appears  that  Che  initial  assumption  that  two  Intervals  would  be 
preferable  to  one  interval  may  be  false.  However,  a closer  examination 
of  the  table  shows  that  AGAUS10S  yielded  substantially  larger  scattering 
fractions  (away  from  0°)  than  did  AGAUS10,  with  the  latter  code's  values 
being  closer  to  the  "REF.  CASE"  values  at  all  choices  of  A. 

About  the  only  conclusions  which  can  be  drawn  at  this  point  are 
that  the  NI  • 1 version  is  by  far  the  most  efficient  when  only  kext  is 
of  interest,  and  that  the  NI  ■ 2 version  is  preferable,  in  spite 
of  its  slower  speed,  if  high  accuracies  are  needed  for  the  scattering 
fractions  at  large  scattering  angles.  From  another  vantage  point, 
however,  the  discrepancies  between  the  NI  « 1 and  NI  ■ 2 scattering 
fractions  at  9 • 90®  and  9 • 180®  are  an  insignificant  fraction  of  the 
9*0  values. 

A few  similar  comparisons  for  Deirmendjian' s cloud  model  C.3 
at  \ • 0. 7 urn  are  given  in  Table  1.2. 5. 7 


1.2. 5. 5 Conclusions : AGAUS10  has  been  demonstrated  to  offer  substantial 
decreases  in  computation  time,  at  selectable  levels  of  accuracy,  over 
either  code  AGAUS9  or  MIE-2.  That  advantage  is  gained  by  letting  the 
program  itself  determine  the  number  of  radii  to  be  used  to  achieve  a 
given  level  of  convergence  for  a specific  "test  quantity",  rather  than 
relying  upon  a user's  guess.  It  also  has  the  additional  property, 
unlike  MIE-2  and  AGAUS9,  that  the  user  can  obtain  some  idea  of  the 
absolute  accuracies  to  be  expected  from  the  results  of  a particular  run. 
There  is  probably  room  for  further  improvements,  of  AGAUS10,  however. 

The  rather  short  existence  of  the  code  (less  than  four  months  in  its 
present  general  form)  has  not  permitted  the  exploration  of  many  of 
its  aspects — such  as  optimizing  the  number  of  radius  intervals  used  with 
the  various  types  of  distributions,  or  of  modifying  it  to  handle 
mixtures  of  aerosol  components  which  have  different  size  distributions 
as  well  as  optical  properties. 
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Table  1.2. 5. 7 Comparison  of  AGAUS10  Results  for  Cloud  Model  C.3 
(X  ■ 0.7  um)  Using  One  and  Two  Size  Intervals, 
Several  Convergence  Levels,  and  Convergence  Tests  on 
Extinction  Coefficient  and  Scattering  Fractions  at 
90*. 


Scattering  Fractions 


DELTA  I NI  INRADI  I Kext  I CPU(a) 


33  3.042  34.81  5.4626E-2  2.2433E-5  2.5350E-4 


65  3.041  | 72.06  ! 5.4583E-2  2.2435E-5  2.3498E-4 


2 130  3.015 


129  3.041 


513  3.021 


5. 3796E-2  2.2180E-5  1.9798E-4 


5. 4585E-2  2.1382E-5  2.3500E-4 


5. 3985E-2  2.0643E-5  2.0382E-4 


Deirmendjian  280  3.021 


5.396E-2 


2.025E-4 


Notes : 


a)  CPU  times  in  seconds,  using  2®  angular  increments. 

b)  Not  comparable  to  other  runs  due  to  use  of  45°  increments. 

c)  Cases  for  which  dashes  ( — ) appear  were  not  run  for  NI-2 


and  convergence  tests  on  both  K and  the  scattered  intensity 


at  0 - 90°. 


65 

3.041 

68.71 

5. 4583E-2 

2.2435E-5 

2. 3498E-4 

f 

66 

3.008 

71.49 

5. 3569E-2 

2.1668E-5 

1. 7792E-4 

65 

3.041 

72.67 

5. 4 58  3D- 2 

2.2435E-5 

2. 3498E-4 
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1.3  MIE  CODE  COMPARISONS 


1.3.1  Introduction:  As  a part  of  the  search  for  ways  of  speeding  up 
Mle-type  calculations  of  aerosol  extinction  and  scattering,  NMSU  has 
conducted  a series  of  studies  of  the  comparative  computational  speeds 
and  results  using  three  separate  computer  codes.  The  three  codes  were: 
MIEGS  - an  ASL  forward-recursion  routine, 

BMIEGS  - an  ASL  backward-recursion  routine, 
and,  DOWN42  - a Radiation  Research  Associates  downward  recursion 
routine. 

In  the  table  accompanying  this  section  of  the  report,  those  three  routines 
will  be  referred  to  as  M,  B and  R,  respectively. 

In  some  cases,  an  "extended"  version^  of  MIEGS  will  be  denoted  as 
M' . All  computations  were  performed  on  WSMR  Univac  1108  computer  system 
B. 


1.3.2  Range  of  Comparisons:  These  studies  used  Mie  size-parameter 
values  (a  - 2ttt/\)  of  1,10,20,40,60,80,100,200,300,500  and  1000,  and 
complex  indices  of  refraction  (n  ■ m-lk)  in  the  forms  k»0,  ® and  m,  and 

m had  the  explicit  values  ra  ■ 1.2,  1.4,  1.6,  1.8  and  5.0. 

The  quantities  which  were  computed  for  comparison  were  the  Mie 

efficiency  factors  Q . and  Q , and  the  average  intensity  I (I.+I^)/2 

J ext  sea  ave  I ^ 

at  scattering  angles  0®,  45*.  90®,  135®  and  180®.  CPU  times  were 

determined  using  a system  supplied  timing  routine  for  each  combination 

of  at,  m and  k. 


1.3.3  Results : The  results  of  the  above  calculations  are  presented 

in  Table  1.3.1  which  shows  the  values  of  at,  m and  k for  which  all  the 

results  (Q  Q , and  I ) agreed  to  better  than  one  (1)  percent, 
^ext*  xsca’  avg  B r 

For  each  ct,  n and  k combination  the  identification  letters  M,  B,  R or  M' 

(see  above)  appear  in  order  of  increasing  computation  time.  If  a 

particular  routine's  outputs  disagreed  with  those  of  the  other  routines 

by  more  than  IX  but  were  within  10X  of  the  other  results,  that  routine 

carries  an  asterisk.  If  a label  (B,  M or  R)  does  not  appear  at  all, 

then  either  it  would  not  run  at  all  (failed  in  a mode  leading  to 

termination  by  the  computer  system  itself),  or  did  not  even  achieve  the 


+The  "extended"  version  of  MIEGS  was  one  in  which  the  dimensions  of  certain 
arrays  were  extended  from  500  to  1500. 
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10Z  agreement  level.  Where  no  entry  at  all  appears  for  some  (a,m,k) 
combinations,  it  means  that  all  three  failed  at  the  machine  detectable 
level.  Finally,  entry  points  showing  only  a single  identifying  letter 
and  a question  mark  represent  cases  in  which  the  other  two  routines 
failed  to  run  at  all  or  gave  such  widely  differing  outputs  that  no 
cross-checking  was  possible. 


1.3.4  Discussion;  In  all  cases  for  which  complete  runs  were  obtainable 
routine  MIEGS  executed  faster  than  BMIEGS  which  in  turn  used  less 
computer  time  than  DOWN42.  BMIEGS  typically  required  twice  as  much 
CPU  time  as  did  MIEGS.  DOWN42  required  from  five  to  twelve  times  as 
much  CPU  time  as  was  needed  by  MIEGS. 

Based  upon  the  results  shown  in  Table  1.3.1,  it  appears  that 
routine  MIEGS  offers  the  best  performance  (accuracy  and  speed)  for  size 
parameters  up  to  about  100  if  the  imaginary  part  k of  the  index  of 
refraction  is  less  than  about  ra/2.  Routine  DOWN42,  however,  appears  to 
be  very  reliable  in  all  the  cases  in  which  it  did  not  fail  at  the 
machine  level,  while  BMIEGS  seemed  to  have  the  most  limited  reliability. 
None  of  the  codes  tested  proved  to  be  capable  of  handling  the  complete 
range  of  parameters  used  in  this  study. 

Since  the  entries  shown  in  Table  1.3.1  for  each  (m,k  and  a) 
combination  are  in  order  of  increasing  computation  times,  that  table 
itself  can  be  used  as  a guide  for  selecting  the  routine  most  likely  to 
be  preferable  for  use  under  various  situations. 


1.3.5  Remarks : As  a part  of  this  study,  NMSU  also  evaluated  the 
so-called  "Deirmendjian  approximations"  [3]  for  the  calculation  of  the 
Q 's.  Those  approximations  were  found  to  be  valid  in  the  regions 
defined  by  Deirmendjian,  but  they  did  not  prove  to  be  faster  in  execu- 
tion than  MIEGS  except  in  the  range  1 < o < 20.  The  computational 
advantages  of  the  Deirmendjian 's  approximations  over  routines  MIEGS 
were  so  small  and  so  restricted  in  applicability  (unusable  for  angular 
intensity  calculations,  for  example)  that  NMSU  sees  little  advantage  to 
their  use  in  any  situation  in  which  MIEGS  can  be  used. 
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1.4  APPROXIMATION  METHODS  FOR  PHASE  FUNCTION  PREDICTION 

1.4.1  Introduction;  In  order  to  surmount  some  of  the  difficulties  in 
determining  the  *true'  or  computed  phase  function  NMSU  has  developed  two 
alternate  methods  of  reducing  the  computational  load.  The  first  method, 
to  be  explained  in  section  1.4.5,  is  essentially  a tabular  one:  given 
a particular  aerosol  model  this  method  will  look  up  precomputed 
scattering  fractions  and  interpolate  over  wavelength  if  necessary.  The 
second  method,  to  be  described  below,  is  the  construction  of  approximate 
phase  functions. 

We  have  used  the  Henyey-Greenstein  phase  function  (hereafter  referred 
to  as  HG) , a new  analytic  phase  function  which  we  shall  refer  to  as 
Goedecke's  phase  function  (hereafter  referred  to  as  GG) , and  a combina- 
tion of  the  HG  and  GG  phase  function  which  we  shall  refer  to  as  the  HG-GG 
phase  function. 


1.4.2  Theory:  The  analytic  HG  phase  function  is  dependent  upon  one 
parameter,  g,  called  the  asymmetry  factor:  g = <cos6>/uj  , where  0 is 
the  scattering  angle  and  H iQ  is  the  albedo  for  single  scattering.  For 
isotropic  scattering  g - 0,  while  for  pure  forward  scattering  g » 1.  In 
our  notation  the  analytic  HG  phase  function  is : 


(1.4-1) 


l wnPn(y),  (1.4-2) 

n*o 


where  P(y)  is  the  approximate  phase  function  at  y = cos0,  and  the 
P (y) * s are  the  Legendre  polynomials.  Using  direct  integration  it  is 
an  easy  matter  to  show  that  g **  5^/(3  u50).  Very  often  the  HG  phase 
function  is  much  too  large  for  angles  close  to  the  forward  or  incident 
directions  (y  # 1)  compared  to  the  computed,  or  'true',  phase  function 
found  by  using  Mie  theory. 

The  GG  phase  function  has  the  following  analytic  form: 


P(y)  - u 


° (1+g2  - 2gy)3/2 


00 

3o  l (2n+l)gn  Pn(y) 
n«o 
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1 
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(1.4-3) 


- Z>0  l (an+l)g"  PR(y) 
n-o 


l 


n-o 


u>'  P (y), 
n n 


(1.4-4) 


where  the  parameters  have  the  same  meaning  as  before;  g^  and  a may  be  chosen 
and  fixed  by  the  following  method.  We  simultaneously  solve  the  following 
equations  for  g^  and  a: 


. - , ot  1 + 81  , ,,  aN  (1-gi) 

P(1)  = P " "o  2 + (1 


V (l-gl)  2f  ’ 


(1.4-5) 


5,  - (L*  - u (a+Dg,  or  2-  - wl  - 2 (2ct+l)g  2,  (1.4-6a,b) 

llo  i Z O l 


where  p = value  of  actual  phase  function  at  y - 1.  In  some  cases 
a > 2.0,  which  means  that  in  these  cases  the  phase  function  can  (and 
sometimes  does)  take  on  negative  values  for  some  y.  This  usually  occurs 
for  y < 0 only,  where  the  phase  function  is  small  compared  to  its  values 
for  y near  1.0. 

The  HG-GG  phase  function  is  a combination  of  the  HG  and  GG  phase 
functions  and  has  the  following  analytic  form: 


U) 


a ^ 

2 (l+g22-2g2y)3/2 


+ (1  - f) 


(l+g22-2g2y)1/2J 


P(y)  - < 


1-82‘ 

u B - 

o . i 


(l+g2  -2g2y) 


3/2 


l^y-p 


y^yS-l. 
(1. 4-7a,b) 


Here,  g0,  a,  3,  y^  are  parameters  that  must  be  fixed  and  are  evaluated 
as  follows:  solve  simultaneously 
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- to  |-  -j 

*(«  5 p Ti^T2  |f(1+g2)  + (1  - f)Cl-82)|.  (1.4-8) 

52  - w0(2a  + l)g22.  (1.4-9) 

Note:  (0^  ■ ajo(cH-l)g2  will  not  work  for  all  cases;  i.e.,  no  real  g2  for 

0 ^ &2  - 1 exists*  The  above  equations  yield  the  same  values  of  g2  and 

a that  we  had  for  g^  and  a for  the  GG  function  alone.  Since  we  expect  that 

for  y<y^  the  actual  phase  function,  and  therefore  the  HG  part  of  the 

HG-GG  phase  function,  is  very  small  compared  to  its  values  for  y>y  , we 
take  these  values  of  g2  and  a as  our  working  values.  We  then  takey..  to 
be  that  value  of  y for  which  the  GG  function  first  goes  negative.  B1is 
found  by  matching  the  GG  part  of  the  phase  function  to  the  HG  phase 
function  at  y - y^.  This  yields 


6 - f + (1  - f)(1  + g2-2gy1)/(l-g2)  > 0.  (1.4-10) 

Therefore  P(y)  > 0 for  all  y and  is  a continuous  function  of  y.  We 

00  — 

must  now  determine  the  to  [P(y)]  - Y to  P (y)  1 : 

n * n 

n-o 

1 

P(y)  P (y) dy.  (1.4-11) 

n 

-1 


to  • 

n 


2n+l 


If  y^>-l,  none  of  the  to^  will  exactly  match  the  computed  or  'true'  <o  . 

Since  it^is  most  important  that  an  approximate  phase  function  have  the 
correct  to  , we  redefine 


to  _ 

P(y)  = =2,  p(y);  then  if  we  write, 

to 

o 

00 

p(y)  - I <on  Pn(u),  we  hav® 

n«o 

1 

wn  - j P (y)  Pn(y)dy  - 
-1 


to 


<0  (^). 
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(1.4-12) 


L . 


<t 


34 
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This  guarantees  that  u)q  ■ &q,  but  uij  will  not  match  oij  exactly.  As 

long  as  the  HG  part  of  the  HG-GG  phase  function  is  small,  the  mismatch 
between  and  ^ will  be  small. 

In  all  the  above  cases  it  is  necessary  to  calculate  only  two  of 
the  coefficients  in  the  Legendre  expansion  of  the  phase  function,  ui 
and  or  uio  and  uij.  This  may  be  done  directly  in  terms  of  the 

an  and  b^  coefficients  of  the  Legendre  series  occurring  in  Mie  theory 

without  evaluating  the  exact  phase  function  of  any  angle. 

We  want  the  phase  function  for  Incident  natural  light  to  satisfy 


P(U)dfl  - 4irui  - 2ir  P(y)dy  - 4ir  -r 
o L 


(1.4-13) 


where  C is  the  total  cross  section  for  scattering  and  extinction, 
respectively.  In  the  notation  of  van  de  Hulst  (1958,  Chs.  2,9), 


So,  if  we  put 


- h | 

2tt 

f 

P(u>  - fy  I F(S,*)d»  . 


P(p)d!l  - K k*  C3C4t, 


and  therefore 


P(y)  - 


F(e,<f>)d<f>. 


ext  o 


Now  F ( 9 , 4>)  ■ i^(8)sini(J)  + 12(6)  cosJ<fc,  hence 


(1.4-14) 


(1.4-15) 


(1.4-16) 


(1.4-17) 


(1.4-18) 
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Si,  I 


I (2n+l) (a  +b  ) 
_ . n n 


where 


This  allows  direct  calculation  of  P(l)  in  terms  of  the  a and  b . We 

now  need  expressions  for  5,  and/or  uL  in  terms  of  the  a and  b . These 

12  n n 

are  given  by  Chu  and  Churchill  [4].  In  terms  of  the  Mie  coefficients 

a , b , these  expressions  are: 


.5  (n(n+!)-3)  2 (2n+l) 

2 n(n+l)(2n+3)(2n-l) 


+ .15  n / * 4.  b*  b ) 

2 (2n+3)KeUn+2an  n+2V 


where  x ■ 2irr/X,  Q 


and  r is  the  particle  radius 
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1.4.3  Comparison  of  Results  for  the  Analytic  Phase  Functions:  Three 
models  have  been  chosen  to  compare  the  analytic  phase  functions  (HG, 

GG  and  HG-GG)  with  the  'true'  or  computed  phase  function.  These  models 
and  selected  wavelengths  at  which  comparisons  were  made  are:  simulated 
white  phosphorus  [1]  at  1.06  and  10.6  microns;  fog  models  WSMRF1,  and 
WSMRF2  at  wavelengths  of  2.5,  10.0  and  11.0  microns:  the  parameters 
for  each  of  the  two  fog  models  are  tabulated  below. 


WSMRF1 

WSMRF2 

Radius,  Minimum 

1.0pm 

0.4pm 

Radium,  Maximum 

40.0pm 

12 . 0pm 

Radius,  Mode 

10. Op 

4.0pm 

Alpha 

3.0 

6.0 

Gamma 

1.0 

1.0 

Particle  Density 

20.0  cm"5 

100.0  cm- 5 

Because  the  GG^  phase  function 

i,  which  uses 

Hi  , w,  and 
o 1 

to  match  the  computed  phase  function,  will  not  always  work,  it  has  been 
abandoned,  and  GG.,  will  be  used  instead  and  referred  to  simply  as  GG:  the 
GG  function  matches  u>q,  ujj  and  P(p-l). 

Values  of  the  analytical  phase  functions  HG  and  HG-GG,  along  with 
the  computed  phase  functions  are  presented  graphically  in  figures  1.4. 3.1 
through  1.4. 3. 8.  A visual  inspection  of  these  figures  shows  that  the 
HG-GG  analytic  phase  function  matches  the  computed  phase  function  far 
better  than  the  HG  phase  function  near  zero  degrees,  and  that  the  HG-GG 
function  apparently  provides  a better  overall  fit  to  the  computed 
phase  function  than  the  HG  function.  Because  the  scales  of  the  graphs 
usually  do  not  permit  sufficient  resolution  near  180#,  the  values  of 
the  computed  and  analytic  phase  functions  are  presented  in  table  1.4. 3.1 
(two  pages)  at  selected  angles.  Table  1.4. 3.1  also  presents  the 
differences  between  the  computed  and  analytic  phase  functions  at  the 
selected  angles,  and  also  the  root  mean  square  error  for  the  points 
considered. 


1.4.4  Discussion:  Inspection  of  table  1.4. 3.1  shows  that  the  HG-GG 
analytic  phase  function  is  superior  to  the  HG  analytic  phase  function 
at  angles  near  zero  degrees.  However,  as  one  approaches  180®  the  HG 
function  sometimes  provides  a better  fit  to  the  computed  phase  function. 
Since  for  most  of  the  distributions  examined  the  backscatter  is  far  less 
than  the  forward  scattering,  it  becomes  important  to  analytically 
fit  the  computed  phase  function  more  closely  near  the  forward  direction 
of  scattering  relative  to  the  backscatter  direction.  Additionally, 
since  the  HG  analytic  function  fits  the  computed  phase  function  best 
near  180®  and  the  HG-GG  analytic  function  is  comprised  of  the  GG 
function  near  0®  and  the  HG  function  near  180®,  we  again  find  that  the 
HG-GG  analytic  function  will  fit  well  near  180®.  The  root  mean  square 
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Figure  1.4. 3.1  Comparison  of  analytic  phase  functions  HG  and  HG-GG  with 
computed  phase  function  for  simulated  white  phosphorous 
at  1.06  pm. 
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Figure  1.4.3. 4 Comparison  of  analytic  phase  functions  HG  and  HG-GG  with 
computed  phase  function  for  fog  model  WSMRF1  at  10.0  van. 
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Figure  1.4. 3.5  Comparison  of  analytic  phase  functions  HG  and  HG-GG  with 
computed  phase  function  for  fog  model  WSMRF1  at  11.0  ym. 
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Figure  1.4. 3. 7 Comparison  of  analytic  phase  functions  HG  and  HG-GG 
with  computed  phase  functions  for  fog  model  WSMRF2 
at  10.0  ym. 
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1.4. 3. 8 Comparison  of  analytic  phase  functions  HG  and  HG-GG  with 
computed  phase  function  for  fog  model  WSMRF2  at  11.0  ym. 
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error  Is  also  presented  for  the  difference  between  the  computed  and 
analytic  phase  functions,  and  again  indicate  thut  the  HG-GG  analytic 
phase  function  provides  a better  fit  to  the  computed  phase  function 
than  the  HG  analytic  phase  function.  Therefore,  it  is  recommended 
that  if  un  analytic  phase  function  is  to  he  used  it  should  be  the  HG-GG 
function.  This  is  particularly  so  if  the  particular  regime  under 
consideration  is  near  the  incident  (forward)  direction. 


1.4.5  Tabular  Method  for  Determination  of  Phase  Function:  Should  the 
user  desire  more  accuracy  than  the  HG-GG  analytic  phase  function  provides, 
it  is  suggested  that  the  following  tabular  method  be  employed.  The 
program  for  creating  such  a table  is  explained  in  Appendix  B. 

In  order  to  reduce  the  computnt tonal  load  in  determining  the 
computed  phase  function  needed  for  the  multiple  scattering  programs.  New 
Mexico  State  University  has  constructed  a tabular  'look  up'  method  for 
certain  fog  models.  The  fog  models  have  been  provided  by  the  Atmospheric 
Sciences  Laboratory,  and  are  based  upon  the  modified  gamma  distribution, 
and  have  been  labeled  WSMRF1  and  WSMRF2  (see  section  1.4.3). 

Program  AGAUSX  has  been  used  to  calculate  values  of  Che  scattering 
fraction  at  2°  Increments  from  zero  to  180°  and  at  selected  wavelengths 
using  the  modified  gamma/ generalized  Khirgian-Mazin  distribution  (type 
5 within  program  AGAUSX).  The  selected  wavelengths  for  both  models  are 
presented  below. 


Selected  Wavelengths  for  Fog  Models  WSMRFl,  WSMRF2 


Wavelength  (pm) 

AX (pm) 

Remarks 

.55 

N.A. 

one  wavelength 

1.06 

N.  A. 

one  wavelength 

2.5 

N.A. 

one  wavelength 

3.0  - 5.0 

.25 

inclusive* 

8.0  - 10. 0 

.25 

inclusive 

10.0  - 12.0 

.5 

inclusive 

e.g.,  3.25,  3.50,  ...  4.75,  5.0 

The  control  program  called  FIND  is  a stand-alone  routine  or  may  be 
incorporated  into  program  ASLSOM  as  a subroutine.  Documentation  for 
program  FIND  may  be  found  in  Appendix  B.  Program  FTND  provides  scatter- 
ing fractions  und  Che  attenuation  coefficient  in  the  proper  format  for 
input  to  ASLSOM  at  the  wuvelengths  given  above  or  will  perform  linear 
interpolation  within  the  wavelength  range  0.55  pm  to  12.0  |im. 

An  error  analysis  has  been  performed  to  determine  the  percent  difference 
between  the  interpolated  phase  function  and  the  computed  phase  function. 
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The  following  table  should  be  used  with  caution,  since  the  percent 
difference  is  at  the  best  only  a rough  indication  of  the  accuracy  of 
the  interpolation  process. 


Percent  Difference  for  Scattering  Fraction  Interpolation 
at  Wavelength  A for  Models  WSMRF1  (E^)  and  WSMRF2 


(e2) 


A(um) 

E X(Z) 

E2(%)  A(um) 

Ixa) 

e2C%) 

1.06 

16 

16  2.5 

40 

58 

3.0 

119 

241  3.25 

100 

123 

3.5 

19 

26  3.75 

22 

9 

4.0 

11 

5 4.25 

14 

1 

4.5 

16 

6 4.75 

16 

5 

5.0 

21 

9 8.0 

8 

7 

8.25 

3 

0 8.50 

5 

2 

8.75 

4 

0 9.0 

5 

1 

9.25 

5 

1 9.5 

6 

2 

9.75 

3 

2 10.0 

6 

4 

10.5 

11 

8 11.0 

30 

29 

11.5 

13 

13 

1.5  ADDITIONAL  STUDIES 

OF  IR  EXTINCTION 

IN  HYGROSCOPIC  AEROSOLS 

1.5.1  Introduction:  In  a 

previous  document 

[1]  some 

results  of  the 

modeling  of 

IR 

extinction 

in  hygroscopic  aerosols  were 

reported.  During 

the  present 

contract  period,  some  additional 

computations  have  been  made 

to  extend  the  ; 

range  of  relative  humidities  beyond  that 

used  previously 

and  to  examine 

the  effects 

of  "hysteresis"  in  the  condensation  and 

evaporation  of  droplets  formed  upon  hygroscopic  condensation  nuclei.  The 

new  calculations  were  carried  out  using  the 

simulated 

white  phosphorous 

smoke  model 

A* 

defined  in 

the  report  referenced  above. 

Computations  were  carried  out  with  an  early  version  of  program 
AGAUS9  at  wavelengths  of  0.55,  1.06,  and  3.6  Um.  Twelve  values  of 
relative  humidity  were  used,  and  mass  accretion  coefficients  were  taken 
from  Hanel  [5]  for  conditions  of  both  increasing  and  decreasing  relative 
humidities. 


1.5.2  Results : The  data  generated  by  the  new  computations  are  presented 
in  graphical  form  in  figures  1.5.1  through  1.5.3.  Figure  1.5.1  shows 
the  absolute  extinction  coefficients  as  a function  of  relative  humidity 
at  X - 0.55  and  X * 1.06  Um.  The  "arrowheads"  indicate  the  direction  in 
which  relative  humidity  is  changing:  An  arrowhead  tilted  toward  the 
right  represents  increasing  relative  humidity,  and  an  arrowhead  tilted 
toward  the  left  represents  decreasing  relative  humidity.  Figure  1.5.1 


Relative  Humidity  (X) 


Figure  1.5.2  Ratio  of  Model  A'  extinction  coefficients  at  A ■ 1.06  ym 
and  A • 3.56  ym  to  values  at  A ■ 0.55  ym  as  a function 
of  relative  humidity.  The  solid  curves  represent 
increasing  relative  humidity,  and  the  dashed  curves 
represent  decreasing  relative  humidity. 
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Figure  1.5 


.3  Attenuation  coefficient  per  unit  mass  of  Wet  Aerosol 
at  0.55  ym  vs  Relative  Humidity  for  NMSU  smoke  model 
A'.  The  solid  curve  is  for  increasing  relative 
humidity  and  the  dashed  curves  represent  decreasing 
relative  humidity. 


clearly  shows  the  "hysteresis"  effect  mentioned  above.  It  will  be  seen 
that  for  relative  humidities  between  about  60%  and  75%,  that  the  value 
of  the  extinction  coefficient  may  vary  by  a factor  of  two  between 
conditions  of  rising  and  falling  relative  humidity.  For  both  wavelengths 
shown,  figure  1.5.1  can  be  interpreted  as  suggesting  that  clearing 
(associated  with  droplet  evaporation)  of  a smoke  cloud  under  conditions 
of  decreasing  relative  humidity  is  slower  than  its  formation  at  any 
given  relative  humidity  between  about  50%  and  75%.  The  figure  also 
demonstrates  that  the  sense  of  changes  in  relative  humidity  may  be  as 
important  as  its  value  in  modeling  hygroscopic  aerosols. 

In  figure  1.5.2  one  will  find  graphs  of  the  ratios  of  the  extinction 
coefficients  at  1.06  ym  and  3.56  urn  to  those  at  0.55  ym  as  a function  of 
relative  humidity,  figure  1.5.2  is  subject  to  various  interpretations. 

One  inference  is  that  the  aerosol  model  used  for  these  calculations  is 
almost  always  (all  relative  humidities  £ 99%)  more  "transparent"  at 
3.56  ym  than  at  0.55  ym,  but  the  same  statement  cannot  be  made  at  1.06  ym. 
This  particular  aerosol,  while  more  transparent  at  10.6  ym  than  at 
0.55  ym  for  relative  humidities  below  about  59%,  is  comparatively 
less  transparent  for  larger  saturation  ratios.  Other  features  discemable 
from  figure  1.5.2  are  that  extinction  at  1.06  ym  as  compared  to  its  value 
at  0.55  ym  is  much  less  sensitive  to  relative  humidity  than  it  is  for 
3.56  ym  radiation.  In  other  words,  the  extinction  at  1.06  ym  can  be 
predicted  more  accurately  from  knowledge  of  the  extinction  at  0.55  ym  than 
can  be  the  extinction  at  3.56  ym  when  no  information  on  relative  humidity 
is  available.  One  further  obvious  inference  from  figure  1.5.2  is  that 
extinction  at  1.06  ym  will  be  greater  than  at  0.55  ym  if  the  smoke  cloud 
is  clearing  because  of  decreasing  relative  humidity. 

Finally,  figure  1.5.3  shows  the  attenuation  coefficient  per  unit 
mass  of  WET  aerosol  at  X - 0.55  ym  as  a function  of  relative  humidity. 

If  that  figure  is  compared  with  figure  1.5.1,  it  will  be  seen  that  the 
attenuation  per  unit  mass  of  the  existing  aerosol  (WET  -*■  nuclei  + accreted 
water)  is  much  less  sensitive  to  relative  humidity  than  the  attenuation 
per  unit  mass  of  dry  aerosol  material.  It  suggests  that  simultaneous 
in  situ  measurements  of  aerosol  mass  concentrations  and  extinction  which 
preserve  accreted  water  should  yield  less  uncertainty  or  variation  in 
extinction  with  changes  in  relative  humidity  than  measurements  made  on 
dry  aerosol  materials  alone. 


1.5.3  Discussion!  The  results  described  above  can,  of  course,  be  applied 
in  detail  only  to  the  aerosol  model  used  in  generating  them  and  at  the 
specific  wavelengths  used.  They  do  suggest,  however,  that  extinction 
can,  under  some  situations  at  least,  vary  by  nearly  a factor  of  two  at 
a given  relative  humidity  between  conditions  of  rising  and  falling 
relative  humidity,  and  that  the  relationship  between  extinction  at  two 
different  wavelengths  may  also  be  dependent  on  the  direction  of  changes 
in  relative  humidity.  Consequently  it  would  seem  that  persons  involved 
in  both  field  measurements  and  numerical  modeling  should  probably  pay  more 
attention  to  meteorological  conditions  than  they  seem  to  have  done  in  the  past. 
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1.6  WAVELENGTH  CHANGE  ALGORITHM 

A feasibility  study  has  been  undertaken  at  NMSU  to  determine 
whether  it  is  possible  to  determine  the  extinction  at  a given  wavelength, 
XG»  from  the  extinction  at  X « .55  ym  for  various  aerosol  distributions. 
To  accomplish  this,  program  AGAUSX  has  been  used  to  find  the  extinction 
coefficient  at  various  wavelengths  for  three  specific  aerosol  models. 

The  models  used  were  WSMRF1,  WSMRF2,  (see  section  1.4.3)  and  a modifica- 
tion of  WSMRF1  which  shall  be  referred  to  as  model  WSMRF3.  Model  WSMRF3 
is  identical  with  model  WSMRF1  with  the  exception  that  for  WSMRF3  the 
mode  radius  was  set  equal  to  8.0  ym  whereas  the  mode  radius  in  WSMRF1 
was  10.0  ym. 

The  results  are  graphically  displayed  in  figure  1.6.1,  where  we 
have  plotted  the  ratio  of  the  extinction  coefficient  at  chosen  wave- 
lengths to  the  extinction  coefficient  at  0.55  ym  against  the  chosen 
wavelength.  The  disparity  between  models  WSMRF1  and  WSMRF3,  which  differ 
only  by  2 ym  in  the  mode  radius  and  are  therefore  related,  highlights 
the  difficulty  in  attempting  to  find  an  analytic  expression  that  would 
be  valid  for  various  aerosol  models  at  sundry  wavelengths.  To  find  an 
analytic  expression  that  could  be  used  for  models  WSMRF1  and  WSMRF2 
(models  that  are  unrelated)  would  be  more  difficult  than  attempting 
to  relate  associated  models. 

Since  aerosol  distributions  may  vary  greatly  in  a short 
period  of  time,  or  may  be  altered  by  battlefield- induced  contaminants, 
it  would  appear  that  estimates  of  extinction  at  various  wavelengths 
derived  from  the  extinction  at  .55  ym  is  risky  at  best.  This  work  is 
also  in  general  agreement  with  the  work  of  Mason  and  Hoidale  [5]  who 
have  assessed  the  utility  of  visibility  as  an  indicator  of  transmittance 
in  the  infrared.  Chylek  [7]  has  investigated  the  dependence  of  the 
extinction  coefficient  against  liquid  water  content,  and  finds  that  there 
is  a definite  wavelength  (X  * 11  ym)  for  which  these  parameters  are 
related.  However  at  different  wavelengths  (e.g.,  0.55  ym)  Chylek  has 
determined  that  there  is  no  unique  relationship  between  the  extinction 
coefficient  and  the  liquid  water  content;  the  extinction  varies  greatly 
as  a function  of  the  size  distribution  when  the  liquid  water  content  is 
held  constant. 
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Figure  1.6.1  Ratio  of  extinction  coefficient  at  various  wave- 
lengths to  the  extinction  at  0.55  Urn  versus 
wavelength  for  fog  models  WSMRF1,  WSMRF2,  and 
WSMRF3. 


Chapter  2 


MULTIPLE  SCATTERING  INVESTIGATIONS 
Introduction 


Two  aspects  of  radiative  transfer  problems  in  which  multiple 
scattering  effects  are  significant  have  been  subjected  to  further 
study  under  this  contract.  The  two  topics  are  (a)  approximation 
methods  aimed  at  speeding  up  calculations  with  concurrent  sacrifices 
in  accuracy,  and  (b)  further  examination  of  the  problem  of  treating 
multiple  scattering  effects  under  illumination  by  beams  with  finite 
cross-sectional  areas.  The  results  of  these  studies  are  presented 
in  more  detail  in  the  following  sections  of  this  report. 


2.1  APPROXIMATION  METHODS 


2.1.1  Introduction:  An  investigation  of  various  methods  for  obtaining 
numerical  solutions  to  the  equation  of  radiative  transfer  by  approxi- 
mate methods  has  been  undertaken  at  NMSU.  There  are  two  basic  approaches, 
the  first  of  which  uses  the  complete  phase  function  but  treats  the 
multiple  scattering  only  approximately;  this  is  useful  for  obtaining 
estimates  of  the  flux  and  albedo.  The  second  method  uses  an  approximate 
phase  function  and  treats  the  multiple  scattering  exactly;  this  is 
useful  when  the  angular  dependence  of  the  intensity  is  desired. 

We  have  reviewed  numerous  methods  for  elongated  phase  functions 
which  occur  when  the  particle  size  is  greater  than  the  wavelength  under 
consideration  that  have  appeared  in  the  literature.  These  include: 
approximating  the  phase  function  by  only  a few  terms  in  the  Legendre 
expansion;  calculating  first  order  scattering  exactly  and  approximating 
higher  order  scattering;  dividing  the  phase  function  into  various 
segments  such  that  an  exact  solution  may  be  found  for  angles  near 
incidence;  truncation  of  the  phase  function;  approximate  analytic  phase 
functions;  similarity  relations  (see  Figure  2.1.1);  and  asymptotic  curve 
fitting.  None  of  the  above  methods  will  work  well  under  all  conditions. 
Additionally  we  have  tried  curve  fitting  to  calculated  intensities 
(see  Figure  2.1.2);  comparing  intensity  distributions  within  a cloud  to 
those  calculated  up  to  a final  optical  depth  that  is  the  same  as  that 
of  the  intensity  within  the  cloud  (see  Figures  2.1.3  and  2.1.4);  and 
interpolation. 
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X Error  v.r.t.  non- isotropic  solution 
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Figure  2.1.3  Comparison  of  internal  vs  'external'  intensities.  The  solid  and  dotted  lines  corres 
pond  to  transmitted  and  reflected  internal  intensities,  respectively.  The  dashed 
and  dash-dot  lines  correspond  to  'external'  intensities.  Where  a comparison  line 
does  not  appear  the  lines  are  coincident.  The  scale  should  be  divided  by  ten  for 
the  scattered  intensities. 
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Figure  2.1.4  Comparison  of  internal  vs  'external'  intensities.  The  solid  and  dotted  lines  corres- 
pond to  transmitted  and  reflected  internal  intensities,  respectively.  The  dashed  and 
dash-dot  lines  correspond  to  'external'  intensities.  Where  a comparison  line  does  not 
appear  the  lines  are  coincident.  The  scale  should  be  divided  by  ten  for  the  scattered 
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While  it  would  appear  that  one  can  always  find  a method  that  will 
adequately  represent  the  angular  distribution  of  intensity  at  any 
optical  depth  and  wavelength,  it  would  be  unusual  if  the  same  method 
will  adequately  describe  the  intensity  distribution  at  other  optical 
depths,  etc.  Therefore,  because  none  of  the  above  methods  work  well 
at  various  optical  depths  and  angles  of  incidence,  we  have  also 
employed  a tabular  method  (see  section  1.4.5). 

In  dealing  with  approximation  methods  for  the  equation  of  radiative 
transfer  we  have  found  that  there  appear  to  be  two  basic  approaches: 

(1)  methods  that  use  the  complete  phase  function  but  only  treat  the 
multiple  scattering  problem  approximately,  and  (2)  methods  that  use  an 
approximate  phase  function  but  treat  the  multiple  scattering  problem 
in  an  exact  manner.  The  first  of  these  is  generally  useful  for 
obtaining  estimates  of  the  albedo  and  flux,  which  contain  no  information 
about  the  angular  dependence  of  the  intensity,  while  the  second  does 
give  the  intensity  as  a function  of  angle,  albeit  only  in  an  approximate 
manner.  As  we  are  primarily  concerned  with  the  angular  distribution  of 
the  intensity,  the  emphasis  in  what  follows  will  be  placed  entirely  on 
the  second  method.  Wherever  possible  a full  set  of  references  will  be 
provided.  It  should  be  borne  in  mind  that  for  any  specific  radiative 
transfer  problem,  an  approximate  solution  may  be  found  that  adequately 
represents  the  original  distribution  of  intensity.  However,  it  does 
not  appear  that  any  one  approximate  method  can  adequately  describe  the 
variation  of  intensity  with  optical  depth,  particularly  for  internal 
intensities,  nor  can  they  accurately  predict  the  variation  of  intensity 
with  other  parameters. 


2.1.2  Methods : 

1)  Exact  solutions  to  the  equations  of  transfer  for 
slightly  elongated  phase  functions. 

Irvine  [8]  has  examined  the  effect  of  approximating  an  elongated 
phase  function  by  one  that  contains  three  terms  or  less  in  the  Legendre 
expansion  of  the  phase  function.  This  was  done  for  situations  with 
azimuthal  symmetry,  conservative  scattering,  and  for  T < 1.  Generally 
he  finds  that  as  the  optical  depth  increases  the  approximated  intensity 
approaches  the  true  value  but  the  process  is  very  slow.  Irvine  finds 
that  the  diffuse  reflection  determined  by  this  approximate  method  "may 
vary  by  a factor  of  2 or  3 from  those  (values)  for  large  particle 
scattering,  while  the  diffuse  transmission  is  qualitatively  different". 

Approximating  an  elongated  phase  function  by  isotropic  scattering, 
where  the  scattering  coefficient  (effectively  the  optical  depth  under 
the  above  conditions)  is  appropriately  reduced  to  some  fraction  of  its 
true  value  (in  order  to  simulate  forward  scattering)  was  also 
investigated  by  Irvine.  Comparison  of  this  method  with  exact  solutions 
showed  that  this  approximation  method  is  inherently  wrong,  producing 
large  errors. 
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2)  Approximate  methods  based  on  the  use  of  exact 
expressions  for  first-order  large  particle  scattering 

Using  the  Neumann  solution  to  the  equation  of  radiative  transfer 
Irvine  [8]  has  calculated  successive  orders  of  scattering  under  the 
conditions  T<1,  w "1*  and  azimuthal  symmetry.  He  then  considers  the 
following:  calculation  of  exact  first-order  scattering  and  approx- 
imating additional  higher  order  scattering  by  simplified  phase 
functions;  i.e.,  isotropic  and  Sobolev.  His  conclusions  based  upon 
this  method  are  as  follows:  the  diffusely  transmitted  intensity  is 
well  represented,  but  the  reflected  intensity  may  be  only  qualitatively 
correct. 

3)  Romanova's  small-angle  approximation 

If  the  phase  function  is  highly  peaked  it  is  possible  to  separate 
that  portion  of  the  intensity  that  changes  most  rapidly  with  angle 
compared  with  the  slowly  varying  component.  Thus,  if  the  phase  function 
is  highly  peaked,  and  the  single  scattering  is  primarily  in  the  forward 
direction,  the  following  approximations  are  made  in  the  equation  of 
transfer.  Writing  the  equation  of  transfer  as 


l. 

- ± J I(T,y')P(y,y')dy't  (2.1-1) 


we  replace  the  factor  y in  the  first  term  of  the  above  equation  on  the 
left-hand  side  by  y and  neglect  backscattering  in  the  boundary 
condition.  Setting  y "1  it  can  be  shown  [9,10]  that  the  exact 


solution  to  this  approximate  problem  is: 

00 

I(T,y)  - { l (2£+l)Pl(y)e 


(2.1-2) 


where  F is  the  flux,  P^(y)  are  the  Legendre  polynomials  and  K^l-uy/ (2H+1)  , 

OO 

where  to^  is  defined  by  the  relation  P(cos9)  - £ u)^  P^(cos0), 

P(cos0)  being  the  phase  function.  Romanova's  method  [11]  is  to  now 
separate  the  intensity  into  two  parts;  that  of  the  small  angle 
approximation  and  that  due  to  the  slowly  varying  part  of  the  intensity 
(called  the  correction  factor  by  Romanova).  Thus  I - I +T;  substituting 

this  into  the  equation  of  transfer  and  using  the  definition  of  the  small 
angle  approximation  we  find  that 


I P dy'  + (y  -y) 

o 


(2.1-3) 
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The  above  equation  may  be  solved  by  any  of  the  numerous  numerical 
methods,  the  primary  advantage  being  the  (supposed)  slow  varying  of  I 
as  a function  of  angle.  This  should  simplify  the  quadrature  problem 
usually  encountered  for  elongated  phase  functions.  Irvine  [12]  has 
examined  the  computations  of  Romanova  [13,14]  and  finds  that  her 
results  agree  with  those  computed  by  Irvine  using  the  doubling  method 
to  within  51  accuracy.  To  the  author's  knowledge  a computational 
code  based  on  the  small  angle  approximation  has  not  been  constructed 
in  this  country,  and  one  must  therefore  rely  upon  Irvine's  feeling 
that  Romanova's  method  "may  effect  a very  considerable  saving  in  time 
over  more  conventional  methods." 


4)  Truncation  of  Phase  Function  and  Approximate  Phase 
Functions 

Hansen  [15]  has  made  comparisons  of  intensities  determined  by 
doubling  with  those  determined  by  approximating  the  phase  function  in 
various  manners.  Specifically  his  approximations  were  made  by  1)  trun- 
cating the  peak  of  anisotropic  phase  functions,  2)  using  the  Henyey- 
Greenstein  phase  function,  and  3)  using  the  Kagiwada-Kalaba  phase 
function,  the  latter  two  being  analytic  functions. 

In  truncating  the  peak  of  anisotropic  phase  functions  Hansen 
has  taken  the  slope  of  the  logarithm  of  the  phase  function  as  being 
constant  for  9 < 20°,  with  the  slope  equal  to  that  of  the  untruncated 
phase  function  at  20®.  Once  the  peak  of  the  phase  function  has  been 
truncated  the  optical  depth  and  co  must  be  appropriately  scaled. 

This  approximation  effectively  replaces  photons  that  are  scattered 
through  a few  degrees  by  photons  scattered  through  0°.  As  would  be 
expected  this  method  produces  errors  of  the  order  of  a few  percent 
in:  a)  angles  near  the  incident  angle,  b)  in  the  direct  backscatter, 
and  c)  in  emergent  angles  near  grazing,  as  photons  moving  nearly 
parallel  to  the  surface  stand  a greater  chance  of  being  scattered  than 
do  photons  that  are  moving  perpendicular  to  the  surface.  Unfortunately 
Hansen  gives  no  graphs  of  the  transmitted  intensity  but  one  can  easily 
see  that  (a)  above  would  be  the  largest  source  of  error  in  the  trans- 
mission. 

In  using  the  analytic  Henyev-Greenstein  (HG)  and  Kagiwada-Kalaba 
(KK)  phase  functions  Hansen  has  found  that  neither  of  the  above  phase 
functions  can  accurately  duplicate  the  angular  distribution  of  intensity; 
the  HG  phase  function  represents  the  intensity  distribution  better  than 
the  KK  function.  The  analytic  phase  function's  representation  of  the 
distribution  of  intensity  becomes  better  the  more  anisotropic  the 
original  phase  function  and  the  thicker  the  atmosphere  becomes. 


5)  Similarity  Relations 

Irvine  [16]  in  an  excellent  review  paper  presents  a method  of 
approximating  the  solution  for  an  anisotropic  phase  function  in  terms 
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of  a known  solution  for  isotropic  scattering.  Such  a transformation 
may  be  valid  if  the  radiation  field  has  been  smoothed,  which  usually 
only  occurs  for  large  optical  depths  or  integration  over  angle.  The 
appropriate  equations  are: 


5*  - (l-g)c3o/(l-g  d^),  (2.1-4) 

To  ■ S„>V  <2-1-5) 

where  the  superscript  * refers  to  isotropic  values.  This  method  was 
tried  and  the  results  may  be  seen  in  Figure  2.1.1.  As  expected  the 
fit  is  not  very  good  and  improves  with  increasing  optical  depth. 


6)  Asymptotic  Curve  Fitting 

van  de  Hulst  [17]  has  given  a method  of  essentially  fitting  a 
curve  to  the  intensity  at  large  optical  depths  and  for  J > 0.6.  His 
method  is  presented  in  a way  that  appears  to  be  easily  programmable. 

In  addition  van  de  Hulst' s method  is  set  up  to  receive  output  of 
reflection  and  transmission  functions  that  have  been  computed  by  the 
doubling  method.  In  reference  16,  van  de  Hulst  implies  that  one  needs 
an  optical  depth  of  about  32  to  apply  his  formulae,  but  in  further 
papers  [18,19]  he  seems  to  imply  that  the  results  of  the  asymptotic 
fitting  are  accurate  for  lower  values  of  the  optical  depth.  Unfortun- 
ately van  de  Hulst  only  provides  formulae  for  normal  incidence,  and 
therefore  independent  of  azimuth. 


2.1.3  Additional  Methods:  A number  of  simple  methods  were  examined 
to  determine  their  feasibility  and  accuracy.  These  consisted  of: 

1)  curve  fitting,  2)  comparison  of  internal  versus  'external' 
intensities,  and  3)  interpolation. 


1)  Curve  Fitting 

As  the  distribution  of  the  intensity  for  single  scattering  follows 
the  'probability  distribution'  of  the  phase  function,  it  was  hoped  that 
one  might  be  able  to  adequately  fit  the  true  intensity  by  using  the 
Henyey-Greenstein  phase  function  normalized  to  some  point  on  the  'true' 
intensity  distribution.  Figure  2.2.2  shows  some  results  of  this  method. 
The  results  for  T«2  are  rather  anomalous  insofar  as  the  fit  is  rather 
good  from  0°  to  36°.  However  the  results  begin  to  deviate  after  this. 
Two  fits  were  tried  for  T ■ .0625,  neither  producing  very  good  results. 
Subsequent  to  this  the  same  procedure  was  computerized  and  fits  were 
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tried  for  various  other  intensity  distributions,  normalizing  at 
sundry  points  on  the  distributions.  The  best  fit  was  obtained  by 
normalizing  the  curves  at  the  maximum  of  the  Intensity  distribution, 
but  the  "average"  error  was  still  between  30  and  100  percent.  Coupling 
the  above  error  with  the  fact  that  apriori  knowledge  of  the  intensity 
distribution  is  needed,  the  method  was  given  up  as  inadequate. 

2)  Comparison  of  internal  versus  'external'  intensities 

Because  the  doubling  method  of  numerically  solving  the  equation 
of  transfer  is  the  fastest  computer  method  a comparison  of  the  results 
from  doubling  were  made  with  those  from  the  Gauss-Seldel  method, 
which  will  produce  internal  intensities.  Thus  the  transmitted  and 
reflected  intensities  at  internal  optical  depths  of  .0625,  0.5  and  a 
final  optical  depth  of  2.0,  as  computed  by  the  Gauss-Seidel  method, 
were  compared  with  the  intensities  produced  by  the  doubling  method  at  final 
optical  depths  of  .0625,  .5  and  2.0.  The  results  presented  for  two  values 
of  to  , the  albedo  for  single  scattering,  may  be  seen  in  Figures  2.2.3 
and  2.2.4.  The  results  are  striking  insofar  as  the  transmitted 
Intensities  from  both  programs  match  exceedingly  well.  This  Is  a 
natural  consequence  of  the  extreme  directivity  of  the  phase  function; 
most  of  the  light  is  scattered  in  the  forward  direction.  Unfortunately 
the  same  cannot  be  said  for  the  reflected  intensities;  agreement  is  only 
achieved  at  the  boundary  (t*0):  as  expected,  the  agreement  is  better  as  u) 
decreases . 


3)  Interpolation 

Interpolation  was  tried  using  the  intensities  determined  by  the 
Gauss-Seidel  method.  This  met  with  little  success  as  there  were  not 
enough  values  of  the  intensity  at  various  optical  depths  to  allow 
accurate,  or  even  semi-accurate,  interpolation.  This  situation  can  be 
changed  by  modification  of  the  Gauss-Seidel  computer  code  to  allow  the 
intensities  to  be  printed  out  at  a few  places  within  the  atmosphere. 


2.1.4  Summary : We  may  summarize  the  previous  methods  by  stating 
that  there  is  apparently  no  approximate  method  that  will  adequately 
describe  the  variation  of  intensity  with  angle  under  all  conditions. 
The  methods  that  have  been  reviewed,  or  tried,  appear  to  be  inadequate 
or  else  require  that  a computer  program  of  some  complexity  be  written. 
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2.2  MULTIPLE  SCATTERING  OF  NARROW  RADIATION  BEAMS  BY  AEROSOLS 
2.2.1  Introduction: 

Steady  state  multiple  scattering  by  typical  homogeneous  aerosols  is 
often  described  by  the  radiative  transfer  (RT)  equation 


9*y(r,s)  + Y I(c,g) 


(4ir) 


p(s,s')I(r,s') , 


(2.2-1) 


where  r ■ (x,y,z)  ■ spatial  coordinates,  s ■ unit  vector,  I(r,s)  - 
total  specific  intensity  at  point  r in  direction  s,  Y = extinction 
coefficient  of  medium,  and  p(s,s')~=  phase  function  of  medium.  The 
specific  intensity  I(j,s)  is  the  radiant  power  at  a given  wavelength 
per  unit  solid  angle  around  s that  crosses  unit  area  parallel  to  s, 
at  point  r.  Eq.  (1)  is  thesimplest  of  a hierarchy  of  equations; 
strictly  speaking,  it  applies  only  to  situations  where  there  is  no 
time  dependence,  no  source  of  radiation  in  the  aerosol,  and  polar- 
ization may  be  neglected.  Many  actual  illuminated  aerosols  satisfy 
these  conditions  to  a good  approximation. 


In  theoretical  calculations,  it  is  usual  to  treat  a model  aerosol 
composed  of  spherical  particles,  because  it  is  only  for  spherical 
particles  that  the  scattering  properties  (y,  p(s,s'))  can  be  calculated 
exactly  from  Mie  theory.  For  a homogeneous  aerosol  containing 
several  components  of  different  optical  type  (particles  having 
different  complex  refractive  indices  m^  ■ n^-lk^)  and  different  size 
distributions  f (a)  and  number  densities  for  each  component,  the 
quantities  (Y,pts,s'))  are  given  by 


- Vn 


i 1 


fi(a)Cext(a)da/ 


f1(a)da  = N Cext, 


p(s,s') 


(4it/y)  X 

i 


f ^ (a) 


C(t)(a, 


s ,s ' )da/ 


f t(a)da , 


where  f (a)da  ■ no.  of  particles  of  type  i with  radius  in  (a,  a+da) , 
(i)  1 

C ..(a)  ■ Mie  theory  extinction  cross-section  of  a particle  of  type  i, 
ext  (1) 

radius  a,  C (a,  s,s’)  ■ Mie  theory  differential  scattering  cross- 
section,  from  direction  s'  to  direction  s,  of  a particle  of  type  i, 
radius  a.  The  last  equality  of  Eq.  (2)  defines  C , the  average 
extinction  cross-section,  where  N ■ £ - total  number  density  of 

particles  in  the  aerosol.  i 


r 

i 
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The  phase  function  is  normalized  so  that  dftg  p(s,s')  ■ 4tt  u)  , where 
(Lo,  the  (average)  single-scattering  albedo,  is  given  by 


u 

o 


Y~l  l N 
i x 


f.(a)  C(i) (a)da/ 
l sea 


,fi 


(a)dr  = C /C  . 

sea  ext 


and  the  last  equality  defines  Cg  , the  average  total  scattering 
cross-section.  Note  well  that  the  extinction  coefficient  y - NC 

ext 

is  the  same  quantity  as  the  Kext  vised  in  other  parts  of  this 
report. 


2.2.2  Incident  Beams  and  Formal  Solutions  of  the  RT  Equation: 

Many  problems  of  interest  involve  one  or  several  aerosol  "clouds"  of 
finite  extent,  acted  on  by  a beam  of  radiation  from  some  source. 

Each  of  these  clouds  may  have  non-uniform  composition;  regions  of 
different  composition  have  different  scattering  properties  (y,  p(s,s')), 
but  each  region  may  be  treated  separately,  and  the  intensities 
matched  at  interfaces.  Therefore,  the  generic  problem  for  calculation 
is  that  of  a single  homogeneous  cloud  of  aerosol,  acted  on  by  a 
specified  incident  beam.  If  the  incident  beam  intensity  is  given  on 
the  cloud  surface,  then  it  may  be  found  inside  the  cloud  as  the 
solution  of 

V~r  Iin(H’?)  + Ylin(~’?)  “ °-  (2.2-2) 

/r\ 

Formally,  if  I^n;(r£,s)  given  at  all  points  rr  on  the  surface  £ of 
the  cloud,  then  the  solution  is 

-y£(r ,s)  (V . 

Iin^~’~^  " e ~ Xin  (£  “ (2.2-3) 

where  £(r,s)  ■ distance  from  r back  along  (-s)  to  the  surface  J (see 
appendix'c."") 

A formal  implicit  solution  of  Eq.  2.2.1  in  the  presence  of  an  incident 
beam  is  then  (see  appendix  C) . 

I(r,s)  - Iln(r,s)  + (4tt)~  1 yj  dE,  e_Ys  jdftg ,p(s  ,s  ’ ) I (r-Cs  ,s ' ) . 

o ~ (2.2-4) 

The  formal  iterative  (Bom  series)  solution  of  this  equation  is 
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I(tl+1)(r,s) 


(An) 


&(r,s) 

d€~e 


-y£ 


d^g  • p (s  » S ' ) 1 


(n) 


(r-Cs,s'),  (2.2-5) 


n » 0,1,  ...  , where  I^°\r,s)  = I (r,s),  and  I(r,s)  * £ I^(r,s). 

(n)  th  n"° 

Here,  i'-  '(r,s)  is  called  the  n order  scattered  intensity.  This 

iteration  solution  is  very  useful  in  calculation  of  first  order 

intensities  for  any  geometry,  but  it  becomes  very  hard  to  handle 

numerically  if  many  orders  of  scattering  must  be  calculated. 


2.2.3  Practical  Methods  of  Solution  of  Transfer  Problems: 

For  situations  involving  a given  aerosol  cloud  and  a given  incident 
beam,  the  computational  problem  is  to  calculate  the  radiant  power 
into  a given  detector;  the  detector  is  characterized  by  its  location, 
orientation,  conical  angle,  and  aperture  area.  Regardless  of  its 
location,  the  detector  is  assumed  to  produce  a negligible  perturbation 
of  the  radiation  field.  Therefore,  the  fundamental  problem  is  to 
calculate  the  total  intensities  I(r,s)  inside  the  homogeneous  cloud, 
given  the  cloud  shape  and  size,  an3  the  incident  beam  intensity  on 
the  cloud  surface. 

In  the  absence  of  any  simplifying  symmetries,  the  most  practical 
accurate  method  of  computation  appears  to  be  the  Monte  Carlo 
technique.  This  technique  is  extremely  time-consuming  if  indeed  a 
map  of  I(r,s)  for  all  (r,s)  is  desired. 

In  the  presence  of  symmetries,  many  other  analytical  techniques 
have  been  used  to  provide  computer  algorithms  for  solving  the  RT 
equation.  The  lower  the  symmetry,  the  slower  the  solution.  Most  of 
these  analytical  techniques  require  at  least  that  the  aerosol  cloud 
be  idealized  to  plane-parallel  translation  invariance.  Within  this 
geometry,  only  the  completely  plane-translation-invariant  problem 
that  results  if  the  incident  beam  is  infinite  in  extent  is  amenable 
to  rapid  digital  computer  solution.  Several  analytic  methods  have 
been  applied  to  this  relatively  simple  problem  (e.g.,  doubling 
(Hansen  [20]),  discrete  ordinate  (Mudgett  and  Richards  [21]; 

Gauss-Seidel  (Herman  and  Brownling  [22])).  The  doubling  method  is 
particularly  rapid,  but  it  yields  only  the  transmitted  and  reflected 
intensities,  not  the  intensities  inside  the  cloud. 

The  finite  incident  beam  (searchlight)  problem  is  of  major  interest. 

Inasmuch  as  modern  guidance  systems  make  use  of  laser  radiation.  Some 
analytic  numerical  techniques  have  been  applied  in  the  special  cases 
of  plane  parallel  geometry  with  isotropic  scattering  (Rybicki,  [23] 

Beckett,  [24]  ) and  with  strongly  forward  peaked  scattering  (Weinman,  [25]; 
1976).  However,  it  is  believed  by  the  authors  of  this  report  that  no 
accurate  computerized  analytic  algorithm  for  the  plane  parallel 
searchlight  problem  with  general  scattering  properties  has  yet  been  developed. 
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During  the  tenure  of  this  contract,  an  attempt  was  made  to  validate 
and  apply  the  doubling  method  computer  program  for  the  searchlight 
problem  that  had  been  developed  under  a previous  contract  [1], 

In  principle  this  program  should  be  applicable  to  the  plane-parallel 
normal  incidence  cylindrically  symmetric  searchlight  problem  for 
an  aerosol  with  arbitrary  scattering  properties.  However,  it  was 
found  that  previous  estimates  of  the  Univac  1108  CPU  time  and  storage 
capacity  that  would  be  required  to  provide  sufficient  accuracy 
were  an  order  of  magnitude  too  low;  a realistic  time  estimate  is 
100  hrs  per  calculation,  rather  than  10  hrs;  storage  requirements 
would  be  over  512K  words,  rather  than  128K.  Any  attempt  to  treat 
beams  at  oblique  incidence,  or  non-plane-parallel  clouds,  would 
Increase  the  computation  time  by  yet  another  order  of  magnitude. 

It  became  necessary  to  conclude  that  the  searchlight  problem  is  one 
of  sufficiently  great  complexity  that  its  accurate  solution  by 
analytic  numerical  techniques  probably  exceeds  the  storage  capacity 
and  practical  time  limitations  of  computers  of  the  Univac  1108  gener- 
ation. Indeed,  this  doubling  algorithm  would  strain  the  capabilities 
of  the  CDC  7600  generation.  Evidently,  even  for  the  searchlight 
problem  of  highest  possible  geometric  symmetry,  (normally  incident 
cylindrically  symmetric  beam  on  a plane-parallel  medium) , Monte  Carlo 
algorithms  are  no  more  time-consuming  than  analytic  ones,  if  high 
accuracy  is  desired. 

However,  it  does  appear  possible,  by  means  of  simple  physical 
arguments,  to  predict  narrow  finite  beam  transmission  and  scattering 
approximately , once  the  corresponding  infinite  beam  results  are 
known.  These  arguments  lead  to  the  conclusion  that,  for  narrow 
Incident  beams  and  small  solid  angle  detectors,  the  Bouguer 
exponentially  attenuated  transmission  is  valid,  and  single  scattering 
is  a good  approximation  in  many  situations,  even  for  optically  thick 
clouds.  Similar  conclusions  have  been  reached  by  several  investi- 
gators (Anderson,  [26]  ; Beckett [24];  Zuev  [27]),  using  methods  different  from 
those  used  here.  In  what  follows,  the  simple  physical  arguments 
and  their  predictions  are  presented. 


2.2.4  Searchlight  Beam  Transmission: 

Consider  the  prototypical  finite  beam  transmission  problem  illustrated  in 
Fig.  2.2.1.  It  is  sufficient  to  consider  this  case  of  highest  possible 
symmetry  (normally  incident  cylindrically  symmetric  searchlight  beam 
on  a plane  parallel  cloud). 

Typically,  the  beam  spread  angle  for  a laser  beam  may  be  taken  as  the 
diffraction  limit,  9 * 0.6X/R^,  where  R.  ■ radius  of  laser  aperture 

< 1mm.  This  ranges  from  6 * .00024  radians  for  visible  light  of 

A ■ 0.4  pm  to  8g  = .006  raalans  for  infrared  radiation  of  A - 10  um. 

The  cloud  thickness  L may  range  from  a few  meters  to  several  kilometers, 
as  may  the  detector  positioning  at  distance  D from  the  incident  cloud  face. 


J 
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Figure  2.2.1  Finite  Beam  Transmission  Geometry. 


The  cloud  optical  depth  yL  mav  vary  from  verv  small  (z  0)  to  very 
large  (10  or  more).  But  multiple  scattering  can  possibly  be 
important  only  for  optically  thick  clouds  (yL  > 0.3).  Consider 
that  an  optical  depth  of  unity  corresponds  to  an  actual  thickness 
L between  10m  (dense  fog)  and  10  km  (light  haze) . Then  y < L 


min 


= 10' lm. 


The  beam  entrance  radius  !U,  area  A0  - tR‘,  and  flux  ^ are  deter- 
mined by  how  far  the  laser  is  situated  behind  the  entrance  face  of 
the  cloud.  (If  the  laser  is  in  the  cloud,  ■ R,,  etc.)  The  beam 
radius  at  any  point  z inside  the  cloud  is  given  by  R(z)  * R + zt?., 
the  beam  area  by  A(z)  - tt(R(z))z,  and  the  beam  flux  at  any  point  z 
by  ■Kz)  - $ A (exp-yz)A(z)  ■ (P  /A(z))exp-y-,  where  P * $ A is  the 
incident  power  in  the  laser  beam1! 


The  detector  is  endowed  wi,th  circular  aperture  A<j,  and  conical 
angle  6.;  typically,  A,  < tt  cm',  and  « 4*  ■ .067  radians.  The 
direct  beam  flux  into  the  detector  is  thus 


$(D)  - * (A  /A(D))exp  - yD.  (2.2-b) 

o o 

Now,  note  that,  in  the  absence  of  the  cloud,  the  direct  beam  flux 
into  the  detector  is 


♦°(D)  - * (A  /A(D)), 
o o 

so  that 

Street  ' <KD)/*°(D)  " exp-yD. 


(2.2-7) 

(2.2-8) 


7Q 


Now  the  multiple  scattering  contribution  (diffuse  flux)  to  the 
detector  must  be  estimated.  It  will  be  necessary  to  calculate 
the  number  of  incident  beam  photons  per  unit  time  that  are  scat- 
tered out  of  the  beam  during  its  travel.  The  number  dn  scattered 
between  z,  z+dz  is  given  by 


dn  = (hv)“  $(z)C  N A(z)dz  ■ (hv)-1  go  A (exp-yz)dz, 

sea  o 03 


(2.2-9) 


where  N ■ no.  of  scatterers  per  unit  volume,  C = (average)  total 
scattering  cross-section  per  particle,  C = C scfc  - single 

— n era  ovt 


scattering  albedo,  C ■ (average)  extinctionCcrols-section  per 


° — - » ext  - — ..  w ^ k'-*- 

particle,  y * NC  t»  ” energy  of  photon.  Therefore,  the  total 

number  of  incident  beam  photons  per  unit  time  that  are  scattered 
out  of  the  beam  across  the  whole  cloud  is 


i 

n = | 


dn  - (hv)-1  go  4>  A (1  - e Y ) . 
o o o 


(2.2-10) 


The  diffuse  power  into  the  detector  is  proportional  to  this  number. 


Now,  imagine  that  an  infinite  incident  beam,  with  flux  $>_,  is  normally 

!/>«*•  AM  A A J J f l_  Am  1 J _ . . _ _ _ •.  t jV  . 1 J 


incident  on  the  same  medium,  with  the  same  detector  positioning  and 
orientation.  Then  it  is  reasonable  to  assert  that  only  those 
incident  beam  photons  that  are  scattered  out  of  a cylindrical  volume 
of  optical  radius  unity  (actual  radius  y-1)  centered  around  the 
symmetry  axis  (where  the  detector  is  situated)  will  contribute 
noticeably  to  the  diffuse  flux  into  the  detector.  Therefore,  the 
number  of  such  photons  per  unit  time,  n , is 


noo  “ (hv)-1  to  $ (ir/y2)(l  - e~YL). 
o o 


(2.2-11) 


The  diffuse  power  into  the  detector  is  roughly  proportional  to  this 
number. 


Furthermore,  suppose  that  the  diffuse  flux  into  the  detector,  ^(D) 
has  been  calculated  accurately  for  this  infinite  beam  case.  Then, 
for  the  corresponding  finite  beam  case,  the  diffuse  flux  into  the 
detector,  4>d(D),  should  be  given  approximately  by 


VD)  - (n/nj  ^(D)  - (yRQ)2  T*lff  , 


(2.2-12) 


where  ^(D)  = Tdi^  $ defines  Tdi^,  the  infinite  beam  diffuse 


transmission  into  the  detector.  Then,  the  total  transmission,  T, 
for  the  finite  beam  case  is  given  approximately  by 


♦(D)  + «.(D)  T*  . * 

T 3 — = exp " TD  * <yro)!  ’ 


T * exp-yD  + TT-1(yR(D))2  T 


diff  * 


(2.2-13) 


(2.2-14) 


According  to  the  ranges  of  numerical  values  determined  above, 

R(D)  - R + D0g  < 10" 3m  + (10 3m)  (6x10" 3)  - 6m,  and  y < 10” 1 m~\  so 
(YR(D))2?tt  < 0.1.  In  most  applications,  more  likely  values  would  be 
D = 100m,  y - 10-2  m“ 1 , 9 = 10-3,  so  R(D)  a 0.1  m,  yR(D)  * 10'3, 

(yR(D))2/7r  < 10~6.  Thus  It  may  be  concluded  that  the  diffuse  part 
of  the  transmission  for  the  finite  beam  case  is  usually  much  less 
than  one-tenth  of  the  diffuse  part  for  the  infinite  beam  case. 

OO 

Since  T^^  is  almost  always  smaller  than  T^irect  f°r  small-angle 
detectors,  up  to  large  optical  depths,  the  diffuse  transmission 
is  therefore  negligible  in  most  cases  of  interest,  and  the  total 
finite  beam  transmission  is  insignificantly  different  from  the 
Bouguer  law  TdirecC  **  exP  “ YD.  even  for  large  D (at  least  yD  = 20). 


2.2.5  Searchlight  Beam  Scattering: 

Consider  the  finite  beam  scattering  problem  illustrated  in  Fig.  2.2.2. 
The  detector  is  located  at  coordinates  (x,o,z),  and  receives  radiation 
at  angles  (4>*0,9). 


Figure  2.2.2  Finite  Beam  Scattering  ueometry. 


72 


The  first-order  scattered  intensity  Iv  ;(x,o,z,0,o)  at  the  detector 
is  fairly  easy  to  calculate  from  equation  (2.2-5),  given  a simple  form 
of  incident  beam.  For  this  first  order  scattering,  the  incident  beam 
spread  is  obviously  unimportant,  so  it  is  convenient  to  choose  the 
following  solution  of  Eq.  (2.2-5), 


iIn<*,y.*.v.,»>  - e<Ro-r ).*«, 


(2.2-15) 


. 2 ■ „2\ 1/2 


where  y - cosQ,  r **  (xz+yz)  , 0 is  the  unit  step  function,  and 

(R  , $ ) are.  the  beam  (radius,  flux)  at  entry.  Then,  from  Eq.  (2.2-5) 

for  tan-1  (x/z)<0<  tt/2, 

I(1)(x,o,z,y,o)  - yj  e~yl  e"Y(z“tu)0[Ro- 1 x-fi.  sin0|]d£. 

o (2.2-16) 

Now,  the  step  function  is  zero  unless  (x-R  )/sin0  < l < (x+R  )/sin0, 
so  O - - O 

(1)  YR  * p(u) 

I (x,o,z,y,o)  * 4tt  sin'9 — exp(-Yz-yx(l-y) /sin0) , 

(2.2-17) 

valid  for  yR  <<1,  which  is  almost  always  true  for  a laser  beam  incident 
on  a physicaf  aerosol. 

Now,  from  Eq.  (2.2-3)  the  corresponding  infinite  incident  beam  is 


I (z,y)  - (2it)-1  $ 6(u-l)exp(-yz) . 
in  o 


Then,  from  Eq.  (5),  the  first  order  intensity  is 


li1}(z,y)  * (4ir)~l  $ p(y)  [exp(-yz)  - exp(-Yz/u)  ]/(l-u)  , 

o 


(2.2-18) 


valid  for  y > 0. 


(2.2-19) 


The  first  order  flux  into  the  detector  is  equal  to  the  first  order 
intensity  times  the  detector  solid  angle,  ir0^,  for  9d<<l>  Therefore, 
the  ratio  of  the  finite  beam  to  the  infinite  beam  first  order  fluxes 
is  given  by  the  ratio  of  Eq.  (2.2-17)  to  Eq.  (2.2-19) 


_ YR0d-y) 
*(1)  Sin0 

Voo 


1-exp (-Yz(l-y)/y)  ’ 


(2.2-20) 
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valid  for  tan-1(x/z)  < 0 < ir/2.  For  estimation  purposes,  consider 
y ■ 0;  then 

» YR  exp(-yx).  (2.2-21) 

oo  o 

Now,  the  ratio  of  the  finite  to  infinite  beam  multiply  scattered  fluxes 
into  the  detector  should  be  similar  to  the  diffuse  flux  ratio  in  the 
transmission  case,  because  these  multiply  scattered  fluxes  are 
determined  primarily  by  the  number  of  scattered  incident  beam  photons 
in  each  case.  However,  the  relevant  scattering  volume  that  contributes 
appreciably  to  the  multiply  scattered  flux  in  the  infinite  beam  case 
is  now  a cylinder  of  radius  y-1  centered  at  the  off-axis  detector 
position.  Therefore,  photons  scattered  from  the  incident  finite  beam 
must  travel  an  extra  distance  x,  on  the  average,  to  get  to  the  detector, 
compared  to  photons  scattered  from  the  corresponding  infinite  incident 
beam.  So  when  the  detector  is  located  a distance  x from  the  finite  beam 
axis,  the  ratio  of  finite  to  infinite  beam  multiply  scattered  fluxes 
should  be  approximately 

^(mult) /^(rnuit)  _ (n/n  )exp_Yx>  (2.2-22) 

°°  00 

where  (n.n^)  are  the  relevant  numbers  of  photons  per  unit  time  scattered 
out  of  the  incident  beams,  given  by  Eqs. (2.2-10)  and  (2.2-11).  So  it  is 
estimated  that  the  total  flux  scattered  from  a finite  beam  into  a 
detector  located  and  oriented  as  in  Fig.  2.2.2  is 

$>(1)  + <$>(mult)*(e~Yx)yRo  ^1)[l+yR0(^mult)/^1))].  (2.2-23) 

-Yx  (1)  (1) 

Since  e ' yR  is  supposed  to  equal  4>  , an  approximate 

expression  for  the  total  diffuse  flux  into  the  detector  in  the  finite 

beam  case  is 

«-  <t(1)  [i+yR0C<I,imult)/^1))],  (2.2-24) 


which  means  that  multiple  scattering  is  less  important  in  the  finite 
beam  case  than  in  the  infinite  beam  case,  as  long  as  yR  <1.  If  yR 
is  very  small,  the  higher  order  scattered  radiation  into  an  off-ax?s 
detector  whose  cone  of  view  contains  a segment  of  the  incident  beam 

may  be  neglected,  unless  the  ratio  of  the  infinite  beam 

fluxes  is  correspondingly  large.  But  this  ratio  is  commonly  of  order 
unity  or  smaller. 


2.2.6  Summary:  The  foregoing  analysis  indicates  that,  for  a narrow  incident 


beam,  multiple  scattering  effects  are  considerably  less  important  than 
for  a very  broad  incident  beam,  in  most  cases.  For  a narrow  incident 
beam,  the  multiple  scattering  correction  to  the  transmission  given  by 
the  Bouguer  law  is  almost  always  negligible.  The  multiple  scattering 
correction  to  the  single  scattering  prediction  of  the  flux  into  an 
off-axis  detector  is  also  negligible  in  many  cases  of  interest.  The 
magnitude  of  the  multiple  scattering  corrections  can  be  estimated 
easily  from  Eqs.  (14) , (24)  once  the  infinite  beam  results  are  known. 
These  latter  results  can  be  obtained  rapidly  by  analytic  numerical 
techniques  for  many  realistic  model  aerosol  clouds. 


2.2.7  Discussion : 

It  is  important  to  point  out  that  the  direct  plus  first  order  diffus 
transmission  into  an  on-axis  detector  aligned  with  a narrow  incident 
beam  represents  undistorted  transmission,  while  (almost)  all  higher 
order  transmission  contributes  to  distortion  of  a modulated  incident 
beam.  This  follows  because  multiply  scattered  photons  travel  a longer 
path.  Similarly,  the  multiply  scattered  flux  into  an  off-axis  detector 
represents  a distortion  of  a signal.  For  example,  suppose  the  incident 
beam  to  be  pulsed,  and  a reflecting  target  to  be  in  the  cloud,  on 
the  beam  axis.  Then,  photons  single  scattered  from  the  aerosol  particles 
into  the  off-axis  detector  arrive  at  the  detector  coincidentally  with 
photons  scattered  by  the  target  only;  but  multiply  scattered  photons 
that  enter  the  medium  at  the  same  time  arrive  later.  Now,  the  detector 
receives  singly  scattered  photons  as  long  as  its  cone  of  view  intersects 
the  beam,  but  it  receives  photons  scattered  just  by  the  target  only  if 
its  cone  of  view  also  intersects  the  target.  Therefore,  for  such  a 
scenario,  for  cases  where  multiply  scattered  photons  contribute  negli- 
gibly, the  signal-to-noise  ratio  (ratio  of  flux  from  target  to  flux  from 
aerosol  single  scattering  into  the  off-axis  detector)  should  always 
increase  noticeably  as  the  detector  orients  toward  the  target,  unless 
the  target  itself  is  a strong  absorber  and/or  a weak  emitter  at  the 
wavelengths  used. 

It  is  also  worth  pointing  out  that  the  approximate  physical  argument 
given  here  can  probably  be  refined  to  compensate  for  idiosvncracies  of 
particular  model  aerosols.  In  particular,  a sharply  forward  peaked 
phase  function  strongly  augments  the  infinite  beam  total  transmission, 
and  reduces  the  scattering,  but  does  not  much  affect  the  arguments 
given  above  that  relate  the  finite  to  the  infinite  beam  fluxes.  Such 
considerations  remain  to  be  studied  in  future  work. 


Chapter  3 


UPDATES  IN  PROGRAM  ATRANX 


Program  ATRANX,  an  NMSU/ASL  atmospheric  transmittance  prediction 
code  developed  under  previous  contracts,  has  been  extended  and  improved 
as  a small  part  of  the  current  research  program. 

1.  The  Van  Vleck-Weisskopf  line  profile  has  been  added  (case 
IPRO  - 5)  to  provide  a line  shape  which  is  more  appropriate  to  milli- 
meter wavelengths  than  those  originally  coded  into  ATRANX. 

2.  The  aerosol  extinction  subroutines  have  been  modified  to 
include  the  changes  in  sizes  of  hygroscopic  aerosol  particles  which 
occurs  with  changes  in  relative  humidity.  The  size  adjustment  formula 
used  is  that  of  J.  Fitzgerald  (Naval  Research  Laboratory) , and  has  the 
form: 


in  which 


[1  + 


0.061 

yU-o 


1/3 


and 


f is  the  fractional  humidity  (0.2  <[  f _<  0.97), 
tq  is  a particle  radius  for  f - 0, 

Y is  Fitzgerald's  "air  mass  characteristic". 


The  quantity  f is  calculated  within  the  programs  from  data  on  water 
vapor  densities.  Hie  parameter,  Y*  however,  must  be  supplied  by  program 
users.  Its  value  is  entered  in  the  position  (and  with  the  same  format) 
that  was  used  for  the  quantity  EMM  in  previous  ATRAN  versions.  EMM  is  no 
longer  an  input  parameter.* 


Because  the  same  accretion  of  liqvid  water  which  causes  size 
increases  also  brings  about  changes  in  indices  of  refraction,  it  has 
also  been  necessary  to  add  information  about  the  optical  constants  of 
liquid  water.  Rather  than  using  additional  input  data  cards,  a new 
subroutine  (WATER)  has  been  added  to  provide  table  look-ups  for  the 
wavelength  range  of  0.2  to  200  pm.  The  optical  parameters  m (real  part 
of  index  of  refraction)  and  k (imaginary  part)  are  adjusted  according  to: 


^ « taken  to  be  1.0  in  the  latest  version  of  ATRAN. 
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and. 


m - mw  + (ma  - \)(r/tQ)  3, 

k - k + (k  - k ) (r/r  )3. 
w a w o 


In  the  two  preceding  expressions,  the  subscript  "w"  represents  pure 
liquid  water,  and  "a"  represents  the  dry  aerosol. 

The  above  changes  require  the  substitution  of  new  routines 
(GASUB4,  NMSME9,  NMSA9 , NMSB9  and  NMSC9)  for  some  older  routines 
(GASUB3,  NMSMIE,  NMSA,  NMSB,  NMSC) , as  well  as  the  addition  of 
subroutine  WATER. 

Suggested  general  choices  for  the  parameter  y (called  GAMMAF 
in  the  program  listings)  are  given  on  "comment  cards"  in  the  source 
deck  for  subroutine  NMSB9. 

3.  Floyd  Herbert’s  [28]  Galatry  profile  codes  (XTV,  XTS,  XKGSER, 
XKASSY,  VOIGTX  and  XKCFRA)  have  been  substituted  for  NMSU  routines 
LINE01  and  LINE02  used  in  generalized  Voigt  profile  runs  (IPRO  - 2). 

4.  The  tape  reading  routine  TTC0R8  has  been  replaced  by  routine 
TTCORE,  which  contains  a new  "compression  option"  ICMPR  * 2.  The  new 
option  permits  line  strength  comparisons  before  deletion  of  molecular 
types. 

5.  Some  portions  of  the  main  program  and  subroutine  HTRN10  have 
been  restructured  to  provide  somewhat  faster  execution.  The  main 
program  (ATRAN10)  has  been  replaced  by  ATRAN14,  and  HTRN10  by  HTRN13. 

6.  The  entire  code  has  been  edited  to  improve  readability  by 
ordering  statement  labels  and  formats  consecutively  within  each  routine. 
Complete  listings  and  a source  deck  for  the  revised  code  will  be 
provided  to  ASL. 
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Appendix  A 

PROGRAMS  AGAUS9  AND  AGAUS10 


A. 1 MIE  THEORY* 

Mie  theory  [29]  predicts  the  scattering  by  and  the  absorption  in  an 
isolated,  discrete,  homogeneous,  isotropic  sphere  of  diameter  D with 
a known  complex  refractive  index  n - m-ik  relative  to  the  surrounding 
medium  and  illuminated  by  monochromatic  radiant  energy  with  wavelength 
X in  the  surrounding  medium.  The  theory  is  given  in  detail  in  standard 
texts  and  need  not  be  repeated  here.  Instead,  those  elements 
of  theory  needed  for  an  understanding  of  the  numerical  algorithms  used 
in  the  ASL  model  are  included. 


Scatterers  attenuate  beams  of  radiant  energy  by  scattering  some  of  the 
energy  into  directions  other  than  the  incident  or  forward  direction 
and  by  absorbing  some  of  the  incident  energy  within  the  body  of  the 
particle.  The  combined  effect  of  pure  scattering  by  the  particle  and 
true  absorption  within  the  particle  is  termed  extinction.  The  amount 
of  extinction,  scattering,  and  absorption  by  a single  particle  is 
given  in  terms  of  corresponding  equivalent  blocking  areas  or  cross 
sections,  C fc,  Cgca,  and  Cgkg,  respectively.  These  cross  sections 
depend  only  on  the  refractive  index  of  the  particle  n ■ m-ik  and  the 
size  parameter  a -2itt/X,  where  r is  the  particle  radius,  and  X is  the 
wavelength. 

The  transmission,  T,  of  a cloud  of  particles  of  geometric  depth,  d,  and 
number  density,  N,  is  given  by 

T - e'T,  (A-l) 

with  the  optical  depth,  T,  given  by 


1 


where 


T - K 
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The  balance  between  loss  by  scattering  and  loss  by  absorption  is 
frequently  characterized  by  the  albedo  of  single  scattering  u>o,  given  by 


This  section  is  partially  based  on  material  taken  from  ECOM  report 
ECOM-5558  by  R.B.  Gomez,  C.  Petracca,  C.  Querfeld  and  G.B.  Hoidale 
March  1975. 
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(A-4) 


A scatterer  with  to  ■ 1 has  no  absorption  and  is  termed  a conser- 
vative scatterer.  °The  albedo  ui  gives  the  probability  that  a photon 
encountering  the  scatterer  will  be  scattered  into  some  direction 
including  the  incident  direction. 

Although  the  extinction  by  a cloud  of  particles  is  correctly  given  by 
Eqs.  (A-l)  and  (A-2) , two  implicit  assumptions  may  lead  to  improper 
use  of  the  equations.  The  optical  depth  T in  Eq.  (A-2)  does  not 
include  losses  caused  by  absorption  in  the  medium  surrounding  the 
particles. 

This  assumption  obviously  breaks  down  at  wavelengths  for  which 
atmospheric  gases  absorb  appreciably.  The  second  assumption  is  that 
scattered  photons  never  return  to  the  incident  direction,  i.e.,  that  there 
is  no  multiple  scattering.  This  effect  becomes  increasingly  important 
as  optical  depths  exceed  T » 0.1. 

A final  caution  should  be  noted  in  regard  to  absorption  within  the 
particle.  Although  absorption  within  the  particle  is  correctly 
determined  by  the  wavelength-dependent  imaginary  part  k of  the 
refractive  index  m,  the  explicit  mechanism  which  causes  the  absorption 
is  usually  not  specified.  Usually  the  absorption  is  joule  heating,  and 
it  is  sometimes  necessary  to  account  for  the  isotropic  black  body  radia- 
tion emitted  by  the  scatterer  when  its  temperature  rises  above  that 
of  its  surroundings.  There  may  also  be  circumstances  when  quantum 
transitions  occur  in  the  scatterer  followed  by  emission  at  or  near  the 
same  wavelengths.  It  is  incumbent  on  the  user  of  the  numerical 
algorithms  presented  here  to  properly  include  these  effects  since  they 
are  not  automatically  accounted  for  in  these  algorithms. 

All  scattering  properties  of  spheres  are  computed  from  m and  k,  and 
through  the  use  of  the  induced  electric  and  magnetic  multipole  moments 
of  the  sphere  a^  and  b^,  respectively.  The  moments  are  given  bv* 


V (net)  41  (ct)-n  yna)  r (a) 

an  4"  (na)  E,  (a)-n  V (not)  £'  (a)  ’ 
n n n n 

and 

n r(nct)  4'  (a)-4»  (not)  V (ct) 

» ^ n n n n 

n n 4*’  (na)  t,  (a)-f  (na)  C («)  * 
n n n n 


The  prime  denotes  differentiation  with  respect  to  the  argument.  The 

4*  (z)  and  £ (z)  functions  are  Ricatti-Bessel  functions  of  the  first 

and  third  k?nd,  respectively,  and  are  related  to  spherical  Bessel 

functions  j (z)  and  n (z)  by 
n n 

"fyfote  that  n is  used  as  a subscript,  an  integer  index  and  a complex 
index  of  refraction  when  it  is  not  a subscript. 
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Vz)  " 2 V2)’ 


and 


where 


and 


£n(z)  - zjn(z)-i  2nR(z)  - ^(z)  + iXn(z) 


jn(z)  " (2z}  ' Jn+l/2(z)’ 


n„(2)  - (jf)U2  »n+1/2(z) 


(A- 7) 

(A-8) 

(A-9) 

(A-10) 


The  function  J +^/7^2^  the  integral  order  Bessel  function; 

the  function  ^+1/2 ^2^  tbe  half  integral  order  Neuman  function. 

The  extinction  cross  section  is  computed  from 


ext 


- 2?  I <2n+1>  Re 


n»l 


n n 


(A-ll) 


and  the  scattering  cross  section  from 


Csca  -5?  £,<2"+1> 
n*l 

The  various  cross  sections  are  the  basic  quantities  used  in  scattering 
problems,  but  they  are  not  the  quantities  usually  computed  directly 
from  Mie  algorithms.  Instead,  it  is  more  convenient  to  compute 
dimensionless  efficiency  factors  Q and  Q , which  depend  on  n,  k, 
and  a,  and  which  are  multiplied  byetne  geoml^rical  sphere  cross  section 
to  obtain  the  true  cross  section  ■ TTriQi.  Thus, 

OO 

««t  ‘ Jl(2tt+V  *e(W-  (A'U) 

and  „ oa 

QSCa  - ^ l (2n+l)[|an|2+|bn|i]  . (A- 14) 

Although  the  cross  sections  account  for  the  energy  removed  from  the 
forward  beam,  they  do  not  give  any  information  about  where  the 
scattered  photons  go.  This  information  is  contained  in  scattering 
amplitudes  and  intensity  factors  which  relate  the  flux  density  scattered 
through  an  angle  0 relative  to  the  incident  flux  density.  There  are 
two  amplitudes,  S^(0)  and  $2(0),  and  intensity  factors  i^(0)  and  i.,(0). 
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LI 


n 

L 


[ 


which  correspond  to  light  respectively  polarized  perpendicular  and 
parallel  to  the  plane  of  scattering  defined  by  the  direction  of 
incidence  and  the  direction  of  scattering. 


The  intensity  factors  are  related  to  the  scattering  amplitudes  by 

(A-15) 


and 


i1(0)  “ ISj^O)  |2  , 

i2(0)  - |S2(0) | 2 . 


(A-16) 


The  amplitudes  come  from  the  multipole  moments  through 


Sl<»  • l ^Ti>  [V„(0>rtnT„<6>'-  <A-17> 

n-1 


and 


S2<«-  I IfSiT  [Vn(e)+Vn(9)I- 

n-1 


(A-18) 


and  angular  factors  it  (9)  and  X (9)  defined  in  terms  of  associated 


Legendre  functions: 


it  (9)  - P*(cos  9)/sin  9; 


(A-19) 


T (0) 

n 


dP^cos  0) 


d0 


(A-20) 


The  amplitudes  have  relative  phase  5 - argS^  - argS-,. 
Alternative  expressions  frequently  used  are 

dP  (cos  9) 


Ve) 


d(cos  9)  ’ 


(A-21) 


and 


d?r  (9) 


TT  (0)  - cos  0* 7T  (9)  - sin  0 • , . , 

n n d(cos  0)  ’ 


where 


P (cos  0)  - — 
n „n 


(cos20-l)n 


(A-22) 

(A-23) 


2 n!  dcos  0 


* I 


These  functions  satisfy  the  following  recurrence  relations: 


TT  (8)  - cos  9 tt  .(0)--Srff  _ (8)  (A-24) 

n ln-1;  n-1  n-1  n-2  ’ 


and 


T (8)  ■ cos8  [ 7T  ( 8)  — tt  (9)  ]-(2n-l)sin28ir  , (8)+t  (8).  (A-25) 

n n n-2  n-1  n-2 


The  scattering  cross  section  measures  the  ability  of  a particle  to 
scatter  light,  and  it  is  to  be  expected  that  Cgca  is  obtained  from 
an  integral  over  the  scattering  intensity  factors.  Equation  (A-12) 
follows  from 


i 

Csca  " { (11(9)  + dcos0‘  (A_26) 

_ l 

Although  the  intensity  factors  themselves  may  be  used  in  scattering 
calculations,  they  are  primarily  suited  for  computing  flux  densities, 
and  it  is  frequently  more  convenient  to  measure  and  compute  scattered 
light  in  terms  of  radiances.  Radiances  do  not  have  the  1/r2  dependence, 
and  it  is  therefore  unnecessary  to  know  the  distance  from  the 
scatterer  to  the  detector  if  the  detector  field  of  view  is  small  and 
is  filled  by  the  scattering  cloud.  The  phase  function  p(0)  gives  a 
radiance  X scattered  into  the  9 direction  in  terms  of  the  radiance 
Iq  incident  on  the  particle. 

The  phase  function  is  dimensionless  and  is  defined  here  as 


p(9) 


2ttC 


ext 


[iL(8)  + i2(8)]  . 


(A-27) 


The  normalized  phase  function  p(8)dfl/4ir  gives  the  probability  of 
a photon  being  scattered  through  an  angle  6 into  an  element  of  solid 
angle  dfJ  ■ d<j>dcos9.  The  integral  of  the  normalized  phase  function 
is  the  single  scattering  albedo  u>  , which  gives  the  probability  that 
the  photon  is  scattered: 


<2 

o 


l 

p(8)d4>dcos0 

l 


4itC 


ext 


t 

[i^CSHijfP)  ]dcos8 


(A-28) 


- 1 


The  phase  function  contains  a sum  over  the  polarization  states 
implicit  in  the  i^  and  i0  intensity  factors,  and  is  thus  unsuitable 
for  describing  the  polarization  of  the  scattered  light. 


The  phase  functions  can  also  be  represented  by  a Legendre 
series : 

n-1 

p(9)  - l 5JlPJl(cos9),  (A-30) 


where  the  Legendre  expansion  coefficients  to.  are  given  by 


1 

S£  * [p(e)  VCOs9)  d<cos9) 

L i 


(A-31) 


and  P^(cos9)  are  the  usual  Legendre  polynomials. 


A. 2 DETAILED  DESCRIPTION  OF  PROGRAM  AGAUS9 


Introduction 


Program  AGAUS9  is  designed  to  calculate  various  scattering  functions 
like  phase  functions  (Mie  theory),  scattering  fractions  (ASLSOM  code), 
extinction  cross-section,  attenuation  coefficient,  etc.,  for  diverse 
natural  and  artifically  created  polydisperse  atmospheric  aerosols. 

The  program  consists  of  subroutines  ANGLE,  GUSET,  AG9PT1,  AG9PT2, 

AG9PT3,  AG9PRT,  MIEGS,  GAUS , VERFY,  and  WATER.  The  organization  and 
operation  of  these  subroutines  is  controlled  by  the  MAIN  program 
AGAUS9. 


A. 2.1  General  Description  of  the  Operation  of  AGAUS9 

The  first  input  card  contains  the  parameters  IDSTP,  NRADI,  IDUMMY, 
NWAVE,  NINDX,  MQRTE,  NCRDS,  NUNIT,  IT,  IANG,  and  ISCAT.  The  various 
modes  of  operations  of  AGAUS9  that  are  possible,  are  controlled  by 
some  of  these  parameters.  If  IANG  ■ 0,  the  subroutine  GUSET  is 
called  to  choose  'IT'  scattering  angles  between  0°  and  180°  for  use 
in  numerical  integrations  using  Gauss-Legendre  (GL)  quadrature.  GUSET 
returns  two  arrays  of  numbers  to  MAIN.  The  array  C (I ) corresponds 
to  the  cosine  of  scattering  angles  and  the  array  H(I)  corresponds  to 
what  are  called  quadrature  weights.  Using  the  array  C(I),  MAIN 
constructs  Legendre  polynomials,  PL  which  will  be  used  later. 

If  IANG-1,  all  the  calculations  of  the  previous  paragraph  are  skipped. 
Instead,  subroutine  ANGLE  is  called.  It  also  returns  arrays  C ( I ) 
and  H(I)  to  MAIN.  In  this  case,  however,  H(I)  corresponds  to  'IT' 
equally  spaced  angles  between  0°  and  180°,  and  C(I)  corresponds  to 
the  cosine  of  these  angles. 

Next,  subroutine  AG9PT1  is  called.  One  or  more  input  parameters  are 
read  at  this  point  depending  on  the  value  of  IDSTP.  AG9PT1  calculates 
and  returns  two  arrays  R(I)  and  F(I)  which  describe  the  normalized 
size  distribution,  and  VOL,  the  average  'dry'  volume  per  particle.  The 
array  R(I)  contains  values  of  'NRADI'  particle  radii,  and  the  array 
F(I)  contains  values  of  the  distribution  for  the  corresponding  R(I). 

For  IDSTP«3,  or  7,8,9,10,11,12  AG9PT1  also  returns  DENS,  the  particle 
number  density  (per  cm3).  To  find  what  is  done  for  other  distribution 
types  see  the  definition  of  DENS  below. 

Before  entering,  in  MAIN,  the  do-loop  Indexed  by  NWAVE,  the  parameters 
WAVE,  DWAVE,  RELHUM,  DENSH,  and  TEMP  are  read.  This  do-loop  allows 
computations  at  several  wavelengths  during  a single  run.  Additional 
looping  options  are  available,  however.  For  example,  by  assigning 
DWAVE  any  value  less  than  10-4,  one  could  use  the  same  do-loop  to 
handle  several  values  of  relative  humidity  (RELHUM)  in  a run. 
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A switch  parameter  LLLL  is  assigned  a value  0 or  1,  in  a manner 
dependent  upon  the  value  of  IDSTP  and  DENSH.  Later,  through  the 
assigned  value  of  LLLL  the  value  of  particle  number  density  is 
chosen  from  three  different  values  for  it  that  may  be  available. 

(See  the  discussion  of  the  parameter  "DENS",  below). 

In  subroutine  AG9PT2  optical  and  physical  data  (index  of  refraction 
and  mass  density)  for  the  aerosol  material  are  read.  For  the  user's 
convenience,  subroutine  WATER  has  been  added  to  provide  the  relevant 
data  for  water  (see  subroutine  WATER) . If  RELHUM  and  growth 
factor  (EMUA)  are  non-zero,  AG9PT2  takes  into  account  the  particle 
growth  and  the  change  in  density  and  indices  of  refraction  of  the 
aerosol  material  due  to  the  absorption  of  water. 

In  AG9PT2,  subroutine  MIEGX  is  called  for  each  'adjusted'  radius. 
MIEGX  returns  to  AG9PT2  the  values  of  extinction,  scattering  and 
back-scattering  efficiency  factors,  and  the  average  intensity  factors 
(I^+iJ/2,  for  'IT'  scattering  angles  chosen  earlier.  These  scattering 
functions  are  weighted  and  integrated  over  the  size  distribution. 

If  the  aerosol  is  a mixture  of  components  having  different  physical 
and/or  optical  properties,  the  above  calculation  is  repeated  for 
each  component  and  various  functions  are  summed  over  NINDX  (the 
number  of  components).  In  the  end  AG9PT2  returns  to  MAIN,  for  each 
wavelength,  the  average  (sum/total  particle  number  density)  extinction, 
scattering,  and  back-scattering  cross-sections  (CTSUM,  CSSUM,  CRSUM) , 
the  average  intensities  (P(J)),  and  the  total  mass  concentration 
(TMASS,  in  gm/cc,  which  includes  the  mass  of  any  liquid  water 
absorbed  by  hygroscopic  aerosols) . 

Next,  subroutine  AG9PT3  is  called.  It  uses  the  numbers  received 
from  AG9PT2  and  returns  to  MAIN  the  total  extinction,  scattering, 
and  back-scattering  coefficients,  scattering  fractions  and  phase 
functions.  AG9PT3  also  calculates  the  albedo  for  single  scattering 
"ALBDO"  ■ £3  » Cara/Cayf.  It  prints  out  all  the  single  wavelength 
results.  Later  In  AG§?T3 , if  IANG  =»  0,  subroutine  GAUS  is  called. 

Using  the  quadrature  weights  H(I)  calculated  earlier,  GAUS 
generates  the  Legendre  expansion  coefficients  <3^  for  the  phase  functions, 
and  places  them  in  array  0L(I) . w^'s  are  then  used  to  reconstruct 
phase  functions,  PC(I)  and  the  routine  computes  the  rms  (root  mean 
square)  deviation  between  the  original  phase  functions  and  the 
reconstructed  phase  functions.  See  the  section  on  subroutine  GAUS 
for  further  details.  If  IANG  - 1,  the  above  calculation  is  skipped. 

If  NWAVE  > 1,  all  the  calculations  of  AG9PT2  and  AG9PT3  are  repeated 
(NWAVE-1)  additional  times.  In  MAIN,  various  scattering  functions 
are  summed  and  divided  by  NWAVE.  Then,  subroutine  AG9PRT  is  called 
to  print  out  all  the  averaged  results.  It  also  calls  GAUS  once  more 
(if  IANG  - 0)  to  generate  the  coefficients  for  the  averaged 
phase  functions. 

There  are  many  other  options  available  which  the  user  will  find 
mentioned  in  this  report.  One  which  is  of  importance  relates  to  control 
parameter  ISCAT.  If  ISCAT  ■ 0,  scattering  fractions  are  printed; 


l 


if  ISCAT  - 1,  they  are  punched  Cor  written  on  NUNIT)  as  well  as 
printed.  If  IANG  - 0,  scattering  fractions  are  neither  printed 
nor  punched. 


SYMBOL 

AMAX 

CATTN 

CRSUM 

CSSUM 

CTSUM 

DENS 

DENSH 

DRYVOL 

DWAVE 

ELWC 

EMM 


Explanation  of  Symbols  used  in  AGAUS9  (MAIN) 
Explanation  or  Definition 


the  largest  Mie  size  found  in  the  aerosol  distribution. 

the  average  (sum/NWAVE)  attenuation  coefficient  in  square 
meters  per  milligram  of  aerosol  material. 

the  back-scattering  cross-section  in  sq.  micrometers  for 
each  wavelength  and  integrated  over  size  distribution. 

AG9PT3  returns  'coefficients'. 

the  scattering  cross-section  in  sq.  micrometers  for  each 
wavelength  and  integrated  over  size  distribution.  AG9PT3 
returns  'coefficients'. 

the  extinction  cross-section  in  sq.  micrometers  for  each 

wavelength  and  integrated  over  size  distribution.  AG9PT3 
returns  coefficients  . 

the  particle  number  density  (number  per  cubic  centimeter) . 
The  value  of  DENS  may  be  supplied  by  the  user  as  an  input 
parameter,  DENSH.  However,  it  will  be  ignored  when 
IDSTP  * 0 or  greater  than  6,  because  these  distribution 
type  carry  pre-determined  values  of  DENS.  In  the  case 
of  IDSTP  » 1,2, 4, 5 the  user  supplied  value  of  number  density 
will  be  ignored  only  if  it  is  less  than  10“ 4 . In  such 
cases  DENS  is  calculated  from  mass  density  and  concentration 
and  average  volume  per  particle,  and  is  represented  by  DENSC 

the  user-supplied  particle  number  density;  units  are 
particles  per  cubic  centimeter.  See  discussion  of  DENS 
(above) . 

the  average  volume  per  particle  of  dry  aerosol  in  cubic 
micrometers. 

the  wavelength  increment  in  micrometers.  See  NWAVE. 

the  liquid  water  content  in  gm/cm3  (used  only  for  cases 
IDSTP-6  and  12). 

the  refractive  index  of  the  surrounding  medium,  taken  to  be 
a vacuum  in  this  program. 


ENWAV 


D FLOAT  (NWAVE). 


AD-A069  449 


UNCLASSIFIED 


NEW  MEXICO  STATE  UNIV  LAS  CRUCES  DEPT  OF  PHYSICS  F/6  4/1 

STUOIES  ON  THE  DEVELOPMENT  OF  ALGORITHMS  FOR  THE  PREDICTION  OF  —ETC(U) 
DEC  78  A MILLER » 6 H 60EDECKE*  R C SHIRKEY  DAAD07-78-C-0063 
NMSU-PHYS-586-78-1  NL 
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SYMBOL 

FSUM 

GNU 

LANG 

IDSTP 

I DUMMY 
IS  CAT 

IT 

KB  ART 

KEXTT 

KSCAT 

LLLL 

LMAX 

MCRTE 

MQRTE 

NCRDS 


Explanation  or  Definition 


the  numerical  result  for  the  integral  over  the  size 
distribution  function  with  respect  to  radius;  it  is  used 
to  normalize  the  distribution  function  to  one  (equivalent) 
particle  per  cubic  centimeter. 

wave  number  in  cm-1. 

■ 0,  for  the  computation  of  phase  functions  at  'IT' 
Gauss-Legendre  quadrature  angles. 

■ 1,  for  the  computation  of  phase  functions  and  scattering 
fractions  at  'IT'  equally  spaced  angles  between  0°  and  180°. 

identifies  the  type  of  aerosol  size  distribution  to  be 
used.  It  can  only  take  values  between  0 and  12.  See 
AG9PT1  for  more  details. 

a spare  input  parameter. 

■ 0,  for  scattering  fractions  to  be  printed.  * 1,  for 
the  scattering  fractions  to  be  printed  as  well  as  punched 
or  written  on  NUNIT  format. 

the  order  of  Legendre  expansion  for  phase  functions  when 
IANG  * 0,  or  the  number  of  equally  spaced  angles  between 
0°  and  180°  when  IANG  - 1. 

the  average  (sum/NWAVE)  back-scattering  coefficient  per 
km,  integrated  over  the  size  distribution. 

the  average  (sum/NWAVE)  extinction  coefficient  per  km, 
integrated  over  the  size  distribution. 

the  average  (sum/NWAVE)  scattering  coefficient  per  km, 
integrated  over  the  size  distribution. 

a switching  parameter  used  to  control  whether  or  not  particle 
number  density  is  to  be  calculated  from  DRYVOL  and  mass 
concentration. 

=3*IFIX  (AMAX) . Integer  estimate  of  the  optimal  order  for 
Gauss-Legendre  quadrature;  used  only  for  diagnostic  message 
print. 

is  not  used  in  AGAUS9  (or  AGAUS10) . 

* 12345,  will  cause  subroutine  AG9PT2  to  print  efficiency 
factors  QT,  QS,  and  QR,  and  the  normalized  distribution 
function  for  every  value  of  radius  used. 

■ 1,  for  punch  (or  write  on  NUNIT)  of  Legendre  expansion 
coefficients  of  phase  functions.  = 2,  for  storing  input  and/or 
output  data,  to  be  used  later  in  a program  attached  to  AGAUS9, 
on  a device  identified  by  NUNIT. 


SYMBOL 


NINDX 

NRADI 

NUNIT 

NWAVE 

OL(I) 

OLT(I) 

OUT(  ) 

PI 

PL(I,J) 

PSUM(I) 

PSUMT(I) 

QATTN 

R(NRADI) 
RELHUM 
SCAT (I) 

SCATT(I) 

TEMP 


Explanation  or  Definition 

the  number  of  aerosol  components  which  will  have  different 
optical  constants,  mass  density  or  mass  concentration. 

the  number  of  values  of  particle  radius  to  be  used  In  the 
calculations.  (In  effect,  points  on  the  radius  vj?  size 
distribution  function  plot.) 

defines  the  device  on  which  output  data  may  be  stored  In 
lieu  of  actual  card  punching;  may  be  used  to  place 
nominal  card  output  into  data  flies  on  tape,  etc.  The 
default  value  (NUNIT-0)  is  4 (card  punch)  for  UNIVAC  and 
7 for  IBM. 

the  number  of  wavelengths  or  relative  humidity  values  to 
be  treated  in  a given  run.  DWAVE  has  to  be  less  than  10~“ 
for  the  latter. 

the  average  Legendre  expansion  coefficients.  OLT(I) /ENWAVE. 

the  total  Legendre  expansion  coefficients  (summed  over  all 
values  of  'WAVE'). 

an  array  used  for  storing  some  of  the  numbers  to  be  printed 
later. 

- 3.1415926535898 

the  Legendre  polynomials  of  order  ( I— 1)  and  argument  C(J).  C(I) 
are  cosines  of  the  scattering  angles. 

the  average  phase  functions  integrated  over  the  size 
distribution. 

the  total  phase  functions  integrated  over  size  distribution 
and  summed  over  wavelength. 

the  attenuation  coefficient  in  sq.  meters  per  milligram  of 
aerosol  material  for  each  wavelength. 

the  radius  of  the  largest  particle  encountered. 

the  relative  humidity  in  percent. 

the  scattering  fractions  as  defined  in  ASL-SOM  for  each 
wavelength.  'IT'  elements  in  the  array. 

the  average  (sum/NWAVE)  scattering  fractions  (ASLSOM) . 

the  temperature  of  the  atmosphere  in  degrees  C. 


SYMBOL 

TMASS 

the 

VOL 

the 

WAVE 

the 

are 

Explanation  or  Definition 


I 
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A. 2. 2 Subroutine  Angle 

Subroutine  ANGLE  is  called  when  IANG  - 1.  It  is  used  to  compute 
the  values  of  'IT*  equally  spaced  scattering  angles  between  0*  and 
180°.  It  places  these  angle  values  in  the  array,  H(I).  It  also 
computes  the  cosines  of  those  angles  and  places  them  in  the 
array,  C(I) . 


U 


SYMBOLS  Explanation 

C(I)  the  array  containing  the  cosines  of  scattering  angles 

in  the  array  H(I). 

H(I)  the  array  containing  'IT'  scattering  angles  (in  degrees). 

IT  is  the  number  of  scattering  angles  chosen  between  0#and  180®. 

PI  - 3. 1415926535898 

RAD  RADS*H(I);  scattering  angles  in  radians. 

RADS  - PI/180. 


L 

li 

L 

I 


.. 


— 
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A.2.3  Subroutine  GUSET 

Subroutine  GUSET  uses  the  Davis  and  Rabinowitz  algorithm  [30] 
to  choose  n values  of  cos0,  (k  * l,2,...n)  between  the  interval 
-1  < cosd  1 and  the  corresponding  values  of  quadrature  weights  a.. 
The  abcissas,  cos0,  are  n zeros  of  Legendre  polynomial  P (cost? 
while  the  weights  a£n  are  given  by  n n 


■ 2(1  ' ' (A-32) 

Initial  estimates  of  the  zeros  are  obtained  from  the  n successive 
zeros  of  the  Bessel  function  (JQ(Jk)  ■ 0)  via 


- cos[jk/«n  +|)a-Kl-(^)l)/4)1/2]. 


(A-33) 


Final  values  of  the  are  found  by  Newton-Raphson  iteration.  The 
tolerance  of  Legendre  polynomial  zeros  is  set  at  10“ 1". 

SYMBOL  Explanation 

AKN(K)  a^n  (Eq.  (A-32)).  ’IT'  elements  in  the  array. 

IT  n;  it  is  the  order  of  Legendre  expansion  for  phase  functions. 

P (N)  Legendre  polynomials,  P (x) . 

n 

PI  3.1415926535898 

X xkn  (Eq. (A-33) ) . 

XKN(K)  cose,  ; P (cose,  ) - 0 within  10“ 1H. 

kn  n kn 

TOL  the  tolerance  of  zeros  of  Legendre  polynomials. 

Z(I)  jk  (Eq . (A-33) ) . I,k  - 1,2,...  IT. 


I 


J 

] 

1 
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A. 2. 4 Subroutine  AG9PT1 

The  main  purpose  of  subroutine  AG9PT1  is  to  generate  two  arrays 
R(J)  and  F(J);  R(J)  contains  'NRADI'  values  of  radius  (particle  size) 
in  micrometers  (ym)  and  F(J)  contains  the  corresponding  values  of  the 
distribution  function.  The  units  of  unnormalized  F(J)  are  cm-5  um- 1 
and  that  of  normalized  F(J)  are  vim-1. 

The  other  quantities  that  are  calculated  are  the  normalization  factor, 
FSUM  and  the  average  volume  per  particle,  VOL  (in  um~3).  AG9PT1 
also  calculates  some  other  quantities  which  vary  from  one  distribution 
type  to  another,  as  indicated  in  the  flowchart. 

Two  parameters  IDSTP  and  NRADI  read  on  the  very  first  input  card 
in  MAIN  are  of  relevance  to  AG9PT1.  IDSTP  identifies  the  type  of 
t size  distribution  function  to  be  used  in  a given  run,  and  NRADI  fixes 

the  number  of  radii  that  will  be  used  to  describe  the  size  distribution 
function.  The  details  of  each  distribution  type  are  given  later.  First 
we  give  the  explanation  of  symbol  which  have  the  same  explanation  for 
all  distribution  types. 


SYMBOL  Explanation  or  Definition 

AVOL  the  average  volume  per  particle  in  cubic  microns  obtained 

via  analytical  integration  over  the  limits  RLOO  and  RHI-». 
That  has  only  been  done  for  IDSTP  ■ 5,6,8,9,10,11,12. 

DENS  the  particle  number  density  in  cm-3. 

DELRD  ■ (RHI-RLO) /RADS ; increment  between  successive  values  of  R. 

DELLR  increment  in  radius  for  the  case  IDSTP  - 0. 

F(J)  the  array  containing  'NRADI'  values  of  the  size  distribution 

function.  See  the  description  of  AG9PT1  for  more  details. 


FSUM 

IDSTP 

IDUMMY 

NCRDS 

NRADI 

NRADII 


the  numerical  integral  over  the  size  distribution  function 
with  respect  to  radius  between  the  limits  RHI  and  RLO; 
used  to  normalize  the  distribution  function. 

identifies  the  type  of  aerosol  size  distribution  to  be 
used.  See  below  for  more  details. 

not  used  in  AGAUS9. 

not  used  in  this  subroutine. 

the  number  of  radius  values  to  be  used  in  describing  the 
size  distribution  function. 

is  same  as  NRADI  except  for  IDSTP  ■ 7,8,9,10,11.  For  IDSTF-7, 
the  input  value  of  NRADI  is  ignored,  for  IDSTP«8,9,10,11  it  is 
doubled. 


SYMBOL 


Explanation  or  Definition 


NU NIT  is  not  used  in  this  subroutine. 

PI  - 3.1415926535898. 

R(J)  the  array  containing  'NRADI*  values  of  radius  (particle 

size)  in  micrometers. 

RADS  - D FLO AT  (NRADI-1). 

RHI  the  largest  value  of  radius  used  (in  urn). 

RLO  the  smallest  value  of  radius  used  (in  urn). 

VOL  the  average  'dry'  volume  per  particle  (in  umJ)  calculated 

numerically. 


IDSTP 

0 


1 


2 


3 


Description  of  Types  of  Distributions 
DESCRIPTION 


This  is  arbitrary  user-supplied  distribution.  'NRADI'  + 1 
cards  have  to  be  read;  the  first  card  contains  RLO  and  DELLR 
(;an).  The  rest  of  the  cards  carrv  the  values  of  F(J)  and  must 
be  in  order  of  increasing  radius  value. 


This  zero-order  log-normal  distribution  and  the  distribution 
function  is  given  by 


F(R)  - 


/2t?  Hog^o)R 


■exp 


$ 


oge(R/R) 


£oge(oV 


1/2' 


R 5 RBAR;  a = SIGMA,  is  the  standard  deviation.  This 
distribution  type_requires  one  input  data  card  to  read 
in  the  values  of  R,  0,  RLO,  and  RHI. 


This  is  called  double  exponential  distribution  and  its 
distribution  function  is  given  by 

F(R)  - QA  exp(-AR)  + (l-Q)B  exp(-BR). 

Q = CUE.  This  distribution  type  requires  one  input  data 
card  to  read  in  the  values  of  RLO,  RHI,  Q,  A,  and  B.  Q 
is  dimensionless  while  A and  B have  units  of  Urn 


This  model  (Deirmendjian's  "Model  C")  does  not  require  any 
input  data  card.  It  carries  fixed  value  of  DENS,  RLO, 
and  DELRD.  RHI  is  determined  by  the  input  parameter  NRADI. 

F(R)  - 450.2  R -*08 

- 2.251*DELRD*R~*  R >.08 
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DESCRIPTION 


The  distribution  function  of  this  model  (Junge  distribution) 
is  given  by 

F(R)  - QR  , 

Q = CUE.  This  distribution  type  requires  one  input  data 
card  to  read  in  RLO,  RHI,  Q,  A. 

The  distribution  function  for  the  Modified  Gamma/Generalized 
Khirgian-Mazin  distribution  is 

y 

F(R)  - Ra  exp[-(|  ) • 2]  , 
c ' 

a = ALF,  y = GAM,  and  Rc  = RC.  One  input  card  is  needed 
to  read  in  RLO,  RHI,  RC,  ALF,  GAM,  and  ELWC.  ELWC  is  not 
needed  for  type  S distribution  and,  therefore  can  be  left 
blank. 


This  size-distribution  model  (NMSU  Fog  or  Cloud  Model) 
is  very  similar  to  type  5,  except  in  that  the  user  must 
supply  one  additional  input  parameter — namely,  the  liquid 
water  content  (ELWC)  in  gm/cc.  This  model  can  be  used  for 
treating  situations  involving  liquid  water  aerosols  like 
clouds  or  fogs.  For  type  6 runs  one  does  not  need 
to  read  in  the  values  of  EMA,  CAYA,  RHOA,  CONC. 

This  distribution  is  essentially  same  as  Junge' s distribution 
(type  4)  except  that  it  has  fixed  parameters.  The  input 
value  of  NRADI  read  in  initially,  is  ignored.  One  input 
card  is  needed  to  read  in  VIS  (visibility  in  km) ; VIS  is 
used  in  calculating  DENS. 

This  is  a fixed  parameter  Continental  Bi-modal  Model.  It  does 
not  require  an  input  data  card  of  type  2 (see  sec.  A. 3. 3). 

This  is  a fixed  parameter  Maritime  Bi-modal  Model.  It  does 
not  require  an  input  data  card  of  type  2 (see  sec.  A. 3. 3). 


This  is  a fixed  parameter  Urban  Bi-modal  Model.  It  does  not 
require  an  input  data  card  of  type  2 (see  sec.  A. 3. 3). 


This  is  a user-supplied  Bi-modal  Model.  This  requires  one 
input  card  to  read  in  FOA,  FOC,  SGA,  SGC,  RBARA,  RBARC. 

Types  8,9,10,  and  11  use  the  sum  of  two  log-normal  distributions 


F(R) 


2 

l 


i-1 


N 


i 


/2i  «.oge(a1)R 


exp 


r L Ao8e(R/yyi 
<_n  ^MaiVJ 


N^  = FOA,  R^  = RBARA,  = SGA  with  similar  meaning  for 
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IDSTP 


DESCRIPTION 


N-,  R.,,  and  a~.  Note  that  in  Type  8,  9,  10  the  values 
or  SGA  and  SCC  are  iog  (a).  For  type  8,9,10,  and  11  the 
total  number  of  radii  used  will  be  twice  as  large  as 
the  input  value  of  NRADI. 

12  This  model  (Marshall-Palmer  Rain  Model)  is  a simple 

exponential  model  which  assumes  an  empirical  relationship 
between  rain  rate  and  droplet  size  distribution  parameters: 

F(D)  - N exp(-AD). 
o 

Nq  - 0.08  cm-*,  and  A - 41  (RN)“°*21  in  which  RN  - RAIN 
is  the  rain  rate  in  mm/hour.  Diameter  D is  in  cm.  The 
corresponding  size  distribution  function  as  a function  of 
radius  R is  given  by 

F(R)  - 2Nq  exp(-2AR). 

This  distribution  requires  one  input  data  card  to  read 
in  RAIN.  The  values  of  RLO  and  RHI  are  fixed  at  0 and  0.5  cm 
respectively.  Due  to  the  limitations  on  the  range  of  Mie- 
sizes  (subroutine  MIEGX)  type  12  usage  is  limited  to  wave- 
lengths of  the  order  of  1 mm  or  larger.  Since  subroutine 
WATER  does  not  contain  optical  data  for  wavelengths  longer 
than  0.2  mm,  type  12  runs  require  the  user  to  supply  the 
values  of  EMA,  CAYA,  and  RHOA  as  if  rain  were  a non-aqueous 
aerosol. 
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Subroutine  AG9PT1  - Simplified  Flowchart 


Receive  IDSTP  and  NRADII  from  MAIN 


Set  NRADI  - NRADII 


Read  and  print  data  for  distribution  type 
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Definition 
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SYMBOL 


CATTN 

the  total  attenuation  coefficient  in  square  meters  per 
milligram  of  'dry'  aerosol  material  for  each  WAVE.  However, 

CTSUM  corresponds  to  'wet'  aerosol. 

CATTNW 

the  total  attenuation  coefficient  in  square  meter  per 
milligram  of  'wet'  aerosol  material. 

CAY 

the  imaginary  part  of  'effective'  refractive  index  of 
'wet'  aerosol.  When  RELHUM  and  growth  factor  are  non-zero 
particles  absorb  water  vapors. 

CAY  - (CAYW  + (CAYA  - CAYW) /A)  . CAY  =*  CAY /EM. 

CAYA 

the  imaginary  part  of  refractive  index  of  dry  aerosol. 

(Assumed  to  be  negative ; do  not  enter  a value  with  a 
negative  sign  on  the  data  card(s)). 

CAYW 

the  imaginary  part  of  refractive  index  of  pure  water  at 
a given  TEMP (DEG  C)  and  WAVE. 

CH 

- FH/(1-FH). 

[ 

CONC 

the  mass  concentration  in  gm/cc  of  a component  of  dry 
aerosol.  It  is  the  number  of  grams  of  dry  aerosol  per 
cubic  centimeter  of  "cloud"  or  "fog",  etc. 

CONCT 

the  total  mass  concentration  in  mg/cc  of  dry  aerosol. 

l 

* . 

) 1 

CRNEW 

■ FACT1*QR;  both  FACT1  and  QR  correspond  to  the  current 
radius  in  the  do- loop  indexed  over  NRADI,  except  when 

L-l;  CRNEW  does  not  exi9t  for  L»l. 

[ 

CROLD 

* FACT1*QR;  both  FACT1  and  QR  correspond  to  the  previous 
radius  in  the  do-loop  Indexed  over  NRADI  except  when 

L*l;  when  L“l,  they  correspond  to  the  current  radius. 

I 

CRSUM 

the  average  (sum/DENST)  back-scattering  cross-section 
in  square  micrometers,  integrated  over  the  size  distri- 
bution, for  each  WAVE. 

1 

CSNEW, 

CSOLD 

1. 

= FACT1*QS . See  CRNEW,  CROLD. 

1 

1 

1 i 

CSSUM 

the  average  (sum/DENST)  scattering  cross-section  in 
square  micrometers,  integrated  over  the  size  distribution, 
for  each  WAVE. 

1 

b 

1 J3 

CTNEW, 

CTOLD 

- FACT1*QT.  See  CRNEW,  CROLD 

j 

■ 

a 

• 
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SYMBOL 

CTNK 

CTSUM 

CTSUMT 

DENS 

DENSC 

DENST 

DQT 

DRYVOL 

DVOL 

ELWC 

EM 

EMA 

EMASS 

EMM 

EMUA 


Definition 


the  extinction  cross-section  in  square  micrometers, 
integrated  over  the  sire  distribution,  for  each  WAVE  and 
each  component  of  the  aerosol. 

the  average  (sum/DENST)  extinction  cross-section  in 
square  micrometers,  integrated  over  the  size  distribution, 
for  each  WAVE. 

the  total  extinction  cross-section  in  square  micrometers, 
integrated  over  the  size  distribution,  for  each  WAVE. 

the  particle  number  density  per  cubic  centimeter.  See 
DENSC  also. 

the  particle  number  density  per  cc  for  each  component  of 
aerosol.  It  is  calculated  from  mass  density  and 
concentration,  and  average  volume  per  particle.  If  LLLL«1, 
the  calculated  value  of  DENSC  is  replaced  by  DENS,  the 
value  of  which  has  been  determined  or  supplied  else- 
where. (Also  see  description  of  the  MAIN  program.) 

the  total  particle  number  density  per  cc.  See  also  DENSC. 

a term  in  the  numerical  integration  over  the  size  distri- 
bution of  the  extinction  cross-section  (CTSUM)  by  the 
trapezoidal  method. 

the  average  volume  per  particle  in  cc. 

a term  in  the  numerical  integration  over  the  size  distribution 
of  th«  volume  of  a sphere  by  the  trapezoidal  method. 

the  liquid  water  content  in  gm/cc. 

the  real  part  of  the  effective  refractive  index  of  'wet* 
aerosol.  See  CAY.  EM  - (EMW  + (EMA  - EMW) / A) . 

the  real  part  of  refractive  index  of  the  dry  aerosol. 

is  'wet'  mass  concentration  in  gm/cc  of  a component  of 
an  aerosol  when  RELHUM  and  growth  factor  are  non-zero. 

EMASS  and  CONC  should  have  same  value  if  growth  factor  is 
zero. 

the  tefractlve  index  of  the  surrounding  medium,  set  equal 
to  1 here  implying  no  surrounding  medium. 

Hanel's  mass  accretion  coefficient  U*  Valuea  must  be 
supplied  by  the  user  and  depend  on  the  type  or  composition 
of  the  aerosol  being  modeled  as  well  as  upon  the  value  of 
the  relative  humidity. 
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SYMBOL 

EMW 

F(I) 

FACT1 

FH 

FQT(I) 

FR 

FVOLD(I) 

IDSTP 

IERROR 

KEXOLD 

KEXT 

KEXTT 

LLLL 


MERROR 

MCRTE 

MQRTE 

NCRDS 

NINDX 


Definition 

the  real  part  of  refractive  index  of  water  at  a given 
TEMR  (DEG  C)  and  WAVE. 

the  normalized  size  distribution  function  (in  micrometers-1) 

- j * F(L)*RT**2,  or  j * F(L-1)*RT1**2. 

the  fractional  relative  humidity  (saturation  ratio). 

N 

- DQT/CTSUM.  In  this  formula  CTSUM  - £ DQT^  and  DQT  - DQT^ 

N - NRADI-2 , NRADI-1,  or  NRADI.  1-1 

- F(L) *DENSC.  DENSC  is  in  cm-s. 

- DVOL/VOL.  See  FQT(I) . 

identifies  the  type  of  aerosol  size  distribution  being 
used  In  the  run. 

a flag  which  is  set  to  unity  if  the  number  of  terms 

reaches  the  maximum  value  allowed  in  the  dimensions  of 

the  Mie-coefficients  a and  b . 

n n 

the  extinction  coefficient  per  km  summed  over  the  number 
of  components  one  less  than  the  current  value  of  the 
running  loop  index  NR. 

the  extinction  coefficient  per  km  for  each  component  of 
the  aerosol. 

the  total  extinction  coefficient  per  km  for  each  WAVE. 

a switch  parameter.  If  LLLL*0,  the  particle  number 
density  (DENSC)  is  calculated  using  mass  density  and 
concentration,  and  the  average  volume  per  particle  (DRYVOL) . 
If  LLLL“1,  a pre-calculated  or  pre-supplied  value  of 
DENS  is  used. 

an  error  counter:  If  MERROR  exceeds  10,  execution 
is  terminated. 

is  not  used  in  this  version  of  the  routine. 

■ 12345,  QT,  QS,  QR  and  FR  are  printed  lor  each  radius. 

not  used  in  this  subroutine,  but  appears  in  a common  block. 

the  number  of  aerosol  components  which  will  have  different 
optical  constants,  mass  density  or  mass  concentration. 
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SYMBOL 

NRADI 

NUNIT 


OLSTAR 

P(J) 

PI 

PSNEW(J) 

PSOLD(J) 

PSUM(J) 

QR 

QS 

QT 

R(L) 

RELHUM 

RFACT 

RHO 

RHOA 

RHOW 

RT 

RT1 

TEMK 

TEMP 


Definition 

the  number  of  points  on  the  radius  vs.  size  distribution 
function  plot. 

defines  the  device  on  which  the  input  and/or  output  data 
may  be  stored  in  lieu  of  actual  card  punching;  may  be 
used  to  place  nominal  card  output  into  data  files  on 
tape  etc.  The  default  value  (NUNIT«0)  is  4 (card  punch) 
for  UNIVAC  and  7 for  IBM. 

not  used  in  AGAUS9. 

an  array  containing  'IT'  average  intensity  factors 
(il+ij) /2  for  each  radius  RT. 

- it  « it  . 

■ P(J)*F(L),  where  P(J)  corresponds  to  LC^  radius;  L _>  2. 

• P(J)*F(L),  where  P(J)  corresponds  to  (L-1)C^  radius 
except  when  L-l;  in  that  case  PSOLD(J)  • P(J)*F(1), 
and  P(J)  corresponds  to  the  Is  radius. 

average  intensity  factors  integrated  over  the  size 
distribution.  They  are  then  summed  over  NINDX  components 
and  divided  by  DENST. 

the  back-scattering  efficiency  factor  for  each  radius, 
the  scattering  efficiency  factor  for  each  radius. 

the  extinction  efficiency  factor  for  each  radius. 

the  array  containing  'NRADI'  values  of  radius  in  micrometers. 

the  relative  humidity  in  percent. 

- (R(L)-R(L-l) ) * DENSC. 

the  specific  density  (in  gm/cc)  of  'wet'  aerosol. 

the  specific  density  (in  gm/cc)  of  dry  aerosol. 

the  specific  density  (in  gm/cc)  of  water  at  the  temperature 
TEMP  (in  DEG  C) . 

the  radius  of  a particle  after  taking  into  account  its 
growth  due  to  absorption  of  water.  RT  ■ AC*R(L)  - BC/AC. 

- AC*R(L-1)  - BC/AC.  See  RT.  For  L-l,  RT1  - RT. 
the  temperature  of  the  surrounding  in  degrees  Kelvin, 
the  temperature  of  the  surrounding  in  degrees  Centigrade. 
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SYMBOL 

TMASS 

TVOL 

WAVE 


Definition 


the  total  mass  concentration  in  gm/cc  of  wet  aerosol. 

Is  used  to  pass  the  value  of  DRYVOL  from  MAIN  to  the 
subroutine. 

the  wavelength  in  micrometers  at  which  all  the  scattering 
functions  are  computed. 


Subroutine  AG9PT2  - Simplified  Flowchart 


RETURN 
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A. 2. 6 Subroutine  WATER 

The  purpose  of  subroutine  WATER  is  to  relieve  the  user  of  programs 
AGAUS9  or  AGAUS10  of  the  task  of  looking  up  and  keypunching  data  on  the 
index  of  refraction  and  mass  density  of  liquid  water.  This  routine 
receives  the  wavelength  (un)  and  temperature  (®K)  from  the  calling 
program  as  variables  WVD  and  TEMPD.  It  returns  the  mass  density 
through  variable  DENSD,  and  the  real  and  imaginary  parts  of  the  refractive 
index  of  liquid  water  through  the  variables  EMD  and  CAYD,  respectively. 

The  data  on  optical  constants  coded  into  routine  WATER  were  taken  from 
the  tabulation  by  Irvine  and  Pollock  [31],  and  the  water  densities 
were  taken  from  a copy  of  the  C.R.C.  Handbook  of  Chemistry  and  Physics. 

The  optical  data  cover  only  the  wavelength  range  0.20  urn  to  200  urn 
and  are  not  entered  at  uniform  wavelength  increments.  Values  of  the 
real  (m)  and  imaginary  (k)  parts  of  the  refractive  index  at  wavelengths 
lying  between  values  for  which  data  are  entered  in  the  tables  are 
estimated  through  straight-line  interpolation.  Linear  interpolation 
is  also  used  between  tabulated  temperatures  in  calculating  the  mass 
density  (py  ■'  DENSD). 

Methods  Used 

Subroutine  WATER  conducts  separate  binary  searches  of  the 
wavelength  table  [LAMBDA(  )]  and  temperature  table  [TEMP(  )]  to 
find  the  indices  L and  L+l  which  bracket  the  received  wavelength 
(WVD)  and  temperature  (TEMPD).  It  then  uses  linear  interpolation  to 
get  estimated  values  of  EMT(m)  CAYT(k)  and  RH0DEN(pw). 

The  interpolation  formula  used  can  be  written  as 


(x-x^)  , 

with  y ■ m,  k or  Pw  and 

x ■ wavelength  or  temperature. 


y(x)  - y(*fc)  + 


yt+ry* 

x*+rx* 


0 

£ 

n 


"bracket"  means. 


for  example: 


LAMBA(L)  < WAVE  < LAMBDA  (L+l) 


i 
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SYMBOL 

CAYD 

CAYT 

DENSO 

EMD 

EMT 

H 

L 

LAMBDA ( ) 

NSUBI(  ) 
NSUBR(  ) 

P 

POINT 
TEMP(  ) 

TEMPO 

TMCHUR 

RHODEN 

WAVE 


Explanation  or  Definition 

imaginary  part  of  refractive  index  k [double  precision  format]. 

the  interpolated  value  of  k in  single  precision  form;  used 
intermediately  to  hold  summed  quantities. 

interpolated  result  for  the  mass  density  of  liquid  water  - p 
[double  precision  form]. 

real  part  of  index  of  refraction  - m [double  precision  form]. 

the  interpolated  result  for  the  real  part  of  the 
refractive  index  (single  precision  form). 

an  indexing  (integer)  parameter. 

an  Integer  indexing  parameter. 

Array  of  wavelengths  at  which  data  are  entered  for  m and  k; 
[typed  as  "real"]. 

array  of  values  for  k (or,  nimag^nary^ • [typed  as  "real"], 
array  of  data  entries  for  m (or,  nrea^)  [typed  as  "real"], 
an  (Integer)  indexing  parameter, 
an  (Integer)  indexing  parameter. 

array  of  temperature  values  (*K)  at  which  entries  for  py 
exist  [single  precision  form]. 

temperature  at  which  p is  to  be  found  [double  precision 
form]. 

temperature  at  which  value  of  p is  desired  [single  precision 
form]. 

Interpolated  result  for  0w  [single  precision  form] 

single  precision  version  of  wavelength  at  which  values 
of  n and  k are  desired. 


i 


*1 

.1 


J 

% 

* 

A 


| 

J 
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A.  2 . 7 Subroutine  M1EGX 

Subroutine  MIEGX  computes  various  efficiency  factors,  and  intensity 
factors  i^  and  i2  for  each  complex  refractive  index  m and  size 
parameter  a.  The  Ricatti-Bessel  functions  and  their  derivatives  in 
Eqs.  (A-5)  and  (A-6)  are  computed  by  the  forward  recursion  method. 

The  initial  values  used  in  forward  recursion  are 

^o(z)  ■ sin  z , 
sin  z 

(z)  cos  z, 

1 z 


^(z)  ■ cos  z,  and 


Xi  “ + sin  z. 

1 z 

[Note:  Cn(z)  - ^(z)  + i^Cz)  ] 

The  angular  functions  and  Tn  are  also  computed  by  forward  recursions 

from  Eqs.  (A-24)  and  (A-25).  The  initial  values  used  are  it  (0)  ■ 0, 

tt,  ( 9)  ■ 1,  T (0)  - 0,  and  T.  ( 0)  - cos0.  0 

10  1 

The  Mie  series  is  terminated  either  when  two  successive  terms  have 
( | Re(a  ) | + | Im(a  )|  + |Re(b  )|  + |lm(b  )|)  < 10~s,  or  when  the  number 
of  terms  exceeds  ?8+Fct) . F ?s  1.2  for  8 < 51  and  is  1+2. 26cT°  • 6 1 3 for 
ot  > 51.  The  tolerance  value  of  10“ s could  be  decreased  for  higher 
precision,  if  the  user  so  desires. 


The  subroutine  computes  and  stores  arrays  of  an  and  bn  until  convergence 
is  reached  and  then  generates  necessary  tt  and  T functions.  Finally 
it  computes  values  of  Q , Qgca,  Q (aBsorption) , i , i2,  Q (back 
scatter),  p(0),  and  radiation  pressure  Q . For  AGAUS9,  only  tne  values 
of  C)  and  p(0)  are  returnlS. 
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SYMBOL 

ALPHA 

C(K) 

CAY 

EM 

EN 

EYE1(K) 

EYE2(K) 

FACT 

FAN(N) 
FBN(N) 
GAMMA 
I ERROR 

ISW1 

MCRTE 

01STAR 

P(K) 

PIE 

PIN 

REAN 

REBN 


Explanation  or  Definition 

Mie  size  parameter,  a » 2irr/\. 

the  array  of  cosine  of  the  scattering  angles.  There  are 
'IT'  elements  in  the  array. 

the  ratio  of  the  Imaginary  part  to  the  real  part  of 
adjusted  refractive  index.  See  AG9PT2. 

the  real  part  of  adjusted  refractive  index. 

Floating  point  representation  of  N 

i.  (Eq.  A-15)  at  angles  « Arc  cos[C(K)].  'IT'  elements 
in  the  array. 

i,  (Eq.  A-16)  at  angles  • Arc  cos[C(K)].  'IT'  elements 
in  the  array. 

determines  the  cutoff  criterion  to  terminate  the  Mie  series. 
It  is  equal  to  1.2  if  a _<  51  and  is  1+2. 26a  #,*ls  for  a > 51. 
See  Mie  theory  text. 

■ Im  («n)  J (A-5). 

* Im(b  );  Eq.  (A-6). 
n 

the  true  Imaginary  part  of  adjusted  refractive  index. 

a flag  which  is  set  to  unity  if  the  number  of  terms  (N) 

reaches  the  maximum  value  allowed  in  the  dimensions  of  a 

and  b . n 

n 

a switch  parameter  used  in  applying  cutoff  criterion, 
is  not  used  in  this  subroutine. 

■ u).  ; first  order  coefficient  for  Legendre  expansion 
of  Che  average  intensity  P(K). 

■ (1^+1t)/2,  the  average  intensity  at  angles.  Arc  cos  (C(K)). 
There  are  'IT'  elements  in  the  array. 

-it  - 3.1415926535898. 

- 1^(0);  Eq.  (A-17) . 

- Re(a  ) ; Eq.  (A-5). 

n 

- Re(b  );  Eq.  (A-6). 

n 
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SYMBOL 

RN 

RNL1 

RPN 

SGA 

SGMAS 

SGMP 

SGR 

SGS 

SGT 

SN 

SNL1 

SPN 

SUMI1I 


SUMI2I 


SUMI1R 


Explanation  or  Definition 

- Re(^n(na));  Eq.  (A-7). 

- Re(ii-o(na));  Eq.  (A-7). 

- Re(ii)'  (na)) . Prime  indicates  differentiation  with 
respec?  to  the  argument. 

the  absorption  efficiency  factor. 

the  average  value  of  Cos(9),  where  9 is  the  scattering 
angle;  not  used  by  AGAUS9. 

the  radiation  pressure  factor  Q^-not  used  by  AGAUS9. 
the  back-scattering  efficiency  factor, 
the  scattering  efficiency  factor, 
the  extinction  efficiency  factor. 

■ Im(i^  (na));  Eq.  (A-7). 

n 

- Im(tji  (na));  Eq.  (A-7). 

o 

■ Im(ii>' (na)) . Prime  indicates  differentiation  with 
respect  to  the  argument. 

the  imaginary  part  of  the  scattering  amplitude  S^(9); 


Im  Sl(9) 


1 [Im(a  )tt  (P)-flm(b  )t  (9)  ] . 


n(n+l) 


n n 


the  Imaginary  part  of  the  scattering  amplitude  S-,(9); 
io  si(9>  ■ 


the  real  part  of  the  scattering  amplitude  S^(9); 

*•  si<a)  • i $58- 


the  real  part  of  the  scattering  amplitude  S2(9); 


SUMI2R 


- I$^R*<VVe>+R*<*n)T„(,”>- 
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SYMBOL 

SUMRI 

SUMRR 

SUMS 

SUMS1 

SUMS  2 

SUMT 

TAUN 

TERMN 

TO 

TOL1 

TPN 

UN 

UNL1 

UPN 


I 


Explanation  or  Definition 

N 

- I (-l)n(2n+l)Im(a  -b  ). 

n n 

n-1 

N 

- y (-l)"(2n+l)Re(a  -b  ). 

n-1  n n 

- £(2n+l)(|a  I1  + |b  !*). 

n n 

- I -^-^^^[ReCa  )Re(a  )+Re(b )Re(b  ) 

, n n n+1  n n+1 

n-2 

+ Im(a  )Im(a  ..)+Ia(b  )Im(b  ..)]. 
n n+i  n n+l 

N n . 

- y [Re (a  )Re(b  ) + Im(a  )Im(b  )]. 

~~  n^n-1;  n n n n 

n*-: 

- y (2n+l)  Re (a  + b ). 

n n 

- T (0);  Eq.  (A-17) . 

n 

- |Re(an)|  + |lm(an)|  + | Re(b  ) | + | ImCb^) | . TERMN 
determines  the  cutoff  criterion  for  terminating  Mie  series. 

- <Pn(a);  Eq.  (A-7) . 

- ^Q(ot) ; Eq.  (A-7). 

- (a) ; prime  indicates  differentiation  with  respect 

to  ?he  argument. 

- x^Ca);  Eq.  (A-8). 

- ^(a);  Eq.  (A-8). 

- \^(a);  prime  indicates  differentiation  with  respect  to 
the  argument. 


", 


.) 

J 

I 


A. 2.8  Subroutine  AG9PT3 


Subroutine  AG9PT3  receives  the  values  of  various  cross-sections  and 
the  integrated  average  intensity  factors  from  AG9PT2  and  converts 
them  into  more  directly  useful  quantities: 

Coefficient  (per  km)  ■ 10“ 3*DENS*cioss-section, 

X2 

Scattering  fractions  ■ -^2*DENS*intensity  factors, 

X2 

and  Phase  function  ■ — *intensity  factors. 

"ext 

X is  in  ym.  It  then  prints  out  all  the  single  wavelength  results. 

It  also  computes  albedo  for  single  scattering  using  Eq.  (A-4) . It 
then  calls  on  subroutine  GAUS  which,  among  other  things  (see  GAUS  for 
details) , also  returns  a value  for  albedo  calculated  differently. 

The  two  values  are  compared,  and  if  they  differ  by  more  than  0,01 
percent  the  user  is  advised  to  rerun  the  program  using  a larger  value 
of  IT. 

SYMBOL 
ALBDO 
C(I) 

CAY 

CAYNG 
CRSUM 

CSSUM 

CTSUM 

DENS 
EM 


Explanation  or  Definition 

■ CSSUM/CTSUM;  albedo  for  single  scattering.  See  Mie  theory. 

the  cosines  of  scattering  angles.  'IT'  elements  in  the 
array. 

the  ratio  of  the  imaginary  part  to  the  real  part  of  adjusted 
refractive  index  of  the  last  component.  See  AG9PT2. 

- - CAY. 

the  average  back-scattering  cross-section  when  it  is 
received  by  AG9PT3  but  it  returns  to  MAIN  the  value  of  the 
average  back-scattering  coefficient  per  km. 

the  average  scattering  cross-section  when  it  is  received 
by  AG9PT3  but  it  returns  to  MAIN  the  value  of  the  average 
scattering  c efficient  per  km. 

the  average  extinction  cross-section  when  it  is  received 
by  AG9PT3  but  it  returns  to  MAIN  the  value  of  the  average 
extinction  coefficient  per  km. 

the  total  particle  number  density  in  cm-3. 

the  real  part  of  the  adjusted  refractive  index  of  the  last 
aerosol  component. 
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SYMBOL 

EMM 

GNU 

H(I) 

IANG 

IDSTP 

IT 

KRSUM 

KSSUM 

KTSUM 

NCRDS 

NINDX 

NUNIT 

OL(l) 

PFACT 


Explanation  or  Definition 

the  refractive  index  of  the  surrounding  medium.  It  is 
set  equal  to  1 in  AGAUS9. 

the  wave  number  in  cm-1. 

the  array  containing  'IT'  scattering  angles  (in  degrees). 

■ 0,  for  the  computation  of  phase  functions  at  'IT' 
Gauss-Legendre  quadrature  angles. 

“1,  for  the  computation  of  phase  functions  and  scattering 
functions  at  'IT'  equally  spaced  angles  between  0°  and  180°. 

identifies  the  type  of  aerosol  size  distribution  to  be 
used.  See  AG9PT1  for  more  details. 

the  order  of  expansion  for  phase  functions  when  IANG-0, 
or  the  number  of  equally  spaced  angles  between  0°  and 
180°  when  IANG-1. 

the  back-scattering  coefficient  per  km  integrated  over  the 
size  distribution,  for  each  WAVE. 

the  scattering  coefficient  per  km  integrated  over  the  size 
distribution,  for  each  WAVE. 

the  extinction  coefficient  per  km  integrated  over  the  size 
distribution,  for  each  WAVE. 

■ 1,  for  punch  of  Legendre  expansion  coefficients  of  phase 
functions. 

■2,  for  storing  input  and/or  output  data,  to  be  used  later 
in  a program  attached  to  AGAUS9,  on  a device  identified 
by  NUNIT. 

the  number  of  aerosol  components  which  will  have  different 
optical  constants,  mass  density  or  mass  concentration. 

defines  the  device  on  which  the  input  and/or  output  data 
may  be  stored  in  lieu  of  actual  card  punching;  may  be  used 
to  place  nominal  card  output  into  data  files  on  tape,  etc. 
The  default  value  (NUNIT«0)  is  4 (card  punch)  for  UNIVAC 
and  7 for  IBM. 

the  first  coefficient  in  the  Legendre  expansion  of  phase 
functions.  When  OL(l)  disagrees  with  ALBEDO  by  more  than 
0.01  percent  program  AGAUS9  should  be  rerun  using  a larger 
value  of  IT. 

- WAVE*WAVE/(PI*CTSUM*EMM*EMM) . When  PSUM(J)  is  multiplied 
by  PFACT  we  get  the  original  phase  functions. 


SYMBOL 


Explanation  or  Definition 


PI 

PSUM(J) 

SCAT(J) 

SFACT 


-it-  3.1415926535898. 

the  average  intensity  factors  Integrated  over  the  size 
distribution  when  received  by  AG9PT3,  but  they  are  the 
original  phase  functions  for  each  WAVE  when  returned  to 
MAIN. 

contains  'IT'  values  of  scattering  fractions  for  each 
WAVE. 

- WAVE*WAVE*DENS*10“‘/4*PI*PI.  When  PSUM(J)  is  multiplied 
by  SFACT,  one  gets  the  scattering  fractions  of  ASLSOM  code. 


WAVE 


the  wavelength  in  microns. 
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Subroutine  AG9PT3  - Simplified  Flowchart 


RETURN 
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A. 2. 9 Subroutine  GAUS 


Subroutine  GAUS  computes  Legendre  expansion  coefficients  given  by 
Eq.  (A-31)  numerically.  To  do  that  the  Integral  In  Eq.  (A-31)  is 
replaced  by  summation  as  follows: 


S,  - f p(9.  n)P.(cos6.  )a 


kn  kn  ’ 


where  9^  and  a^  have  the  meaning  as  given  In  subroutine  GUSET. 


Using  the  values  of  coefficients  and  Eq.  (A- 30)  phase  functions 
are  reconstructed.  We  call  them  reconstructed  phase  functions, 

Pc(9  ).  GAUS  then  computes  the  rms  (root  mean  square)  deviation 
between  the  original  phase  functions  p(9  ) and  the  reconstructed 

phase  functions  p (9.  ) as  each  successive  term  is  added  to  the  series 
in  Eq.  (A- 30) . Final?y  it  prints  out  the  values  of  coefficients 
and  rms  deviations. 


SYMBOL 


Explanation  or  Definition 


is  the  array  containing  the  original  phase  function  P(®kn) 


is  the  order  of  Legendre  expansion  for  phase  functions. 


NCRDS 


■ 1,  for  punch  of  Legendre  expansion  coefficients  of  phase 
functions 

■2,  for  storing  input  and/or  output  data,  to  be  used  later 
in  a program  attached  to  AGAUS9,  on  a device  identified 
by  NUNIT. 


NUNIT 


defines  the  device  on  which  the  input  and/or  output  data 
may  be  stored  in  lieu  of  actual  card  punching;  may  be  used 
to  place  nominal  card  output  into  data  files  on  tape,  etc. 
The  default  value  (NUNIT"0)  is  4 for  UNIVAC,  and  7 for  IBM. 


OL(LL) 


the  array  containing  coefficients  <2«. 


PC(I) 


the  array  containing  the  reconstructed  phase  functions 


Pc«W* 


RMS (J) 


the  rms  deviation  between  the  original  phase  functions 
and  the  reconstructed  phase  functions. 


A 

' \ 
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A. 2. 10  Subroutine  AG9PRT 

Subroutine  AG9PRT  prints  out  all  the  averaged  (sum/NWAVE)  results 
essentially  in  the  same  way  as  AG9PT3  does  for  each  wavelength. 

If  IANG*0,  subroutine  GAUS  is  called  to  generate  Legendre  expansion 
Coefficients  which  are  used  to  reconstruct  the  averaged  phase 
functions.  GAUS  also  computes  the  rms  deviation  between  the  original 
phase  function  and  the  reconstructed  phase  functions.  See  the 
section  on  subroutine  GAUS  for  more  details.  If  IANG-1,  the  above 
calculation  is  skipped. 

If  ISCAT-0,  the  averaged  scattering  fractions  are  printed;  if 
ISCAT-1,  they  are  punched  (or  written  on  NUNIT)  as  well  as  printed. 

If  IANG-0,  scattering  are  neither  printed  nor  punched. 


SYMBOL 

C(I) 

CATTN 

H(I) 

IANG 


ISCAT 


IT 


Explanation  or  Definition 


the  cosines  of  scattering  angles.  'IT'  elements  in  the  array. 

the  average  (sum/NWAVE)  attenuation  coefficient  in  square 
meters  per  milligram  of  aerosol  material. 

the  array  containing  'IT'  scattering  angles  (in  degrees). 

- 0,  for  the  computation  of  phase  functions  at  'IT' 
Gaus8-Legendre  quadrature  angles. 

■ 1,  for  the  computation  of  phase  functions  and  scattering 
fractions  at  'IT'  equally  spaced  angles  between  0°  and  180°. 

■ 0,  for  scattering  fractions  to  be  printed. 

■ 1,  for  scattering  fractions  to  be  printed  as  well  as 
punched  or  written  on  NUNIT. 

the  order  of  Legendre  expansion  for  phase  functions  when 
IANG«0,  or  the  number  of  equally  spaced  angles  between 
0®  and  180®  when  IANG*1. 


u 

i: 


KBAKT 


KEXTT 


KSCAT 


NUNIT 


the  average  (sum/NWAVE)  back-scattering  coefficient  per 
km,  integrated  over  the  size  distribution. 

the  average  (sum/NWAVE)  extinction  coefficient  per  km. 
Integrated  over  the  size  distribution. 

the  average  (sum/NWAVE)  scattering  coefficient  per  km. 
Integrated  over  the  cize  distribution. 

defines  the  device  on  which  the  input  and/or  output  data 
may  be  stored  in  lieu  of  actual  card  punching;  may  be  used 
to  place  nominal  card  qutput  into  data  files  on  tape,  etc. 
The  default  value  (NUNIT*0)  is  4 (card  punch)  for  UNIVAC, 
and  7 for  IBM. 


SYMBOL 

NWAVE 

PSUM(J) 

SCATT(J) 


Explanation  or  Definition 

the  number  of  wavelengths  or  relative  humidity  values 
to  be  treated  in  a given  run. 

the  array  containing  the  values  of  average  (sum/NWAVE) 
phase  functions  integrated  over  the  size  distribution. 

the  array  containing  the  values  of  average  (sum/NWAVE) 
scattering  fractions. 
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A. 3 PROGRAM  AGAUS10 

This  appendix  contains  a presentation  of  the  details  of  the 
mathematical  procedures  used  for  the  "halving"  procedures  of  AGAUS10, 
explanations  of  program  symbols  found  In  the  routines  unique  to 
AGAUS10,  and  a summary  of  the  ways  In  which  input  data  for  AGAUS10 
differ  from  those  of  AGAUS9. 

Readers  of  the  following  descriptions  of  the  mathematical  proce- 
dures should  bear  In  mind  that  the  general  AGAUS  single  scattering  code 
computes  a number  of  quantities  associated  with  the  Mle  theory  of  the 
interaction  of  a single  spherical  scatterer  and  then  averages  those 
quantities  using  the  aerosol  size  distribution  function,  f(r),  as  a 
weighting  factor.  The  following  definitions  may  also  be  helpful: 

i i 

The  terms  "component"  or  "aerosol  component"  apply  to  cases  in 
which  the  overall  aerosol  may  be  a mixture  of  different  materials 
whose  optical  properties,  mass  density,  growth  factors,  or  mass 
concentrations  are  not  identical.  The  number  of  "components"  to  be 
used  in  a run  is  specified  by  input  parameter  NINDX.  An  example  of  a 
multi-component  aerosol  might  be  one  which  is  a mixture  of  hygroscopic 
and  non-hygroscopic  particles. 

The  terms  "size-interval"  and/or  "original  interval"  refer  to 

subdivisions  of  the  total  range  of  aerosol  radii  (R  . . • RLO  to 

minimum 

^maximum  ■ RHI)  within  which  the  halving  and  convergence  testing 

procedures  are  applied  independently. 


A. 3.1  The  "Halving"  Method 

Most  averages  that  must  be  done  are  of  the  form 
RHI 

(1)  G = J dr  f (r)  G(r)/FSUM, 

RLO 

where  G(r)  is  some  Mie  theory  quantity  for  a particle  of 
radius  r,  and  a given  optical  type  (refractive  index  m - ra-lk), 
and 


FSUM 


RHI 

r 

dr  f(r).  The  size  distributions  f(r),  RLO<r<RHI, 

RLO 


may  be  given  in  analytic  form,  or  numerically  on  a set 
(R(J)}  of  values  of  r.  They  may  or  may  not  be  normalized 


J 
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to  unity  (FSUM  may  equal  1.0,  or  not).  The  relevant  Mie 
theory  quantities  include  C^xt(r),  Csca(r^«  P(r,u)  for  each 

Gauss-Legendre  u,  etc.  Other  averages  that  must  be  done 
have  somewhat  different  form,  e.g., 

(2)  ^ — | dr  f(r)co  (r)Cext(r)/FSlTN, 

ext  ' 

for  the  coefficients  ui^  of  the  phase  function  expression 

CO 

p(U)  • T id  P (U),  in  terms  of  the  coefficients  (r> 

“ n n n 

n-o 

CO 

of  the  phase  function  p(r,u)  “ V oi  (r)  P (u). 

“ n n 
n-o 

The  above  integrals  over  r must  be  done  numerically.  In  order 
to  minimize  computation  time,  the  number  of  points  in  the  numerical 
integration  should  be  chosen  as  small  as  possible,  consistent  with 
adequate  accuracy.  Each  value  of  r that  is  used  requires  the  full  Mie 
calculation  for  the  quantities  cext(r)*  C (r) , ^or  eac^  U, 

etc.,  these  calculations  absorb  a' large  part  of  the  total  computation 
time.  It  has  been  found  in  this  work  that,  instead  of  fixing  the 
number  of  values  of  r to  be,  say,  500,  adequate  accuracy  results,  in 
many  cases,  for  as  few  as  50  to  100  values  of  r;  this  reduction  allows 
a reduction  in  total  computation  time  by  roughly  a factor  of  4 to  8. 

The  numerical  integration  method  which  was  developed  makes  use  of 
successive  halving  of  intervals  until  a preset  convergence  criterion 
is  met.  This  halving  method  allows  an  initial  set  of  unequal  intervals 
to  be  chosen.  For  example,  suppose  that  the  distribution  function  f (r^ 
is  sharply  peaked  around  some  value  of  r,  and  also  has  a long  tail,  as 
in  Fig.  A. 3. 1. 

Then  it  makes  sense  to  choose  unequal  initial  intervals, 
perhaps  as  shown.  The  numerical  integration  then  proceeds  by  halving 
each  of  the  original  intervals,  until  the  convergence  criterion  is  met. 
The  convergence  criterion  that  was  used  is  the  requirement  that 

|G^n+^  - G^I/Ig^I  < A,  where  A is  some  small  preset  number.  Here, 

G^n^  is  the  value  of  the  Integral  G after  the  nth  halving,  separately 
in  each  initial  interval.  The  trapezoidal  rule  was  used  throughout. 

The  value  of  the  Integrand  at  any  given  point  r is  calculated  Just  once; 
it  does  not  have  to  be  recalculated  on  that  point  after  halving. 
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Figure  A. 3.1  Sketch  of  a hypothetical  size-distribution 


The  logical  flow  through  the  halving  integration  code  is  as 
follows: 

1)  Given  RLO,  RHI,  and  f(r),  RLO  < r < RHI,  in  analytic  form; 
or  f(r)  numerically  on  the  equally  spaced  set  r - R(J), 

J - 1,  NRADI. 

2)  For  the  given  f(r),  RLO,  RHI,  choose  the  number  and  size 
of  the  original  intervals.  That  is,  choose  RR(I),  1-1, 

NLAST  - 2** MIN  + 1,  with  RR(1)  - RLO,  RR(NLAST)  - RHI. 

Note  that  NI  - no.  of  initial  intervals  - 2*MIN  - 1,2,4,8,16 

Then  set  NHALV  • no.  of  halvings  in  each  interval  ■ MAX-MIN, 
and  NK  - maximum  no.  of  points  in  each  interval  - 2**NHALV. 

3)  The  maximum  number  of  points  allowed  is  NMAX  - 1+2**MAX. 

In  this  work,  MAX  - 9 was  used,  so  NMAX  - 513.  Calculate 
R(KK),  F(KK)  - f (R(KK)) , KK  - 1,  NMAX,  on  all  the  points 
R(KK)  which  might  be  used  if  the  halving  goes  all  the  way  to 
NMAX  points.  For  f(r)  given  numerically,  this  is  done  by 
linear  interpolation. 


t 

1 
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RR(I+1) 

4)  Calculate  Cext(I)  - j dr  f(r)Cext(r)  in  eac^  Interval 

RR(I) 


no.  I,  I«l,  NI,  by  halving  of  that  interval  and  the 
trapezoidal  rule,  until  the  convergence  criterion  is  met  in 
that  interval,  or  until  the  halving  has  progressed  all  the 
way  to  NK  points,  the  maximum  number  allowed  in  any  interval. 

RR(I+1) 

At  the  same  time,  calculate  FSUMG(I)  • J dr  f(r),  and 

RR(I) 

all  other  needed  averages 
RR(I+1) 

G(I)  - | dr  f(r)G(r),  on  the  same  set  of  points  used  for 
RR(I) 


C ..(I).  (This  means  that  the  convergence  criterion  is 
ext  NI 

applied  only  to  Cext(D).  Then  CflXt  - V Cext(I), 


NI 

FSUM  - l FSUM(I),  etc.;  and  then  Cext  - C^/FSUMG, 
G • G/FSUMG.  This  last  step  is  equivalent  to  making 


RHI 

dr  f (r) 


RLO 


1.0. 
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A. 3.2  Glossary  for  Subroutines  AGXPT1  and  AGXPT2 


SYMBOL 

RLO 
RHI 
DELLR 
NRADI 
RP 

FF(J) 

MIN 
MAX 
NI 

NLAST 
RR(I) 
NMAX 

I NHALV 

NKG 

RBAR 
SIGMA 

CUE 

A 

B 


Explanation 


min.  value  of  particle  radius  In  any  sire  distribution 
max.  value  of  particle  radius  in  any  sire  distribution 
increment  of  radius  in  arbitrary  (user-supplied)  distribution 
no.  of  radii  at  which  arbitrary  distribution  is  given 
a value  of  particle  radius 

value  of  arbitrary  distribution  function  at  radius 
RP  - RLO  + (J-1)*DELLR,  J-l,  NRADI 
integer 
integer  _>  MIN 

2**MIN  ■ number  of  initial  intervals  in  halving  method  of 
integration 

NI+1  - number  of  radii  RR(I)  defining  the  NI  basic  size 
intervals 

lower  radius  of  I initial  interval.  1-1 , NLAST 
RR(1)  - RLO,  RR (NLAST)  - RHI 

2**MAX+1  ■ maximum  number  of  points  which  may  be  used  in 
halving  method 

MAX-MIN  ■ maximum  number  of  interval  halvings  allowed  for 
each  initial  interval 

2**NHALV  ■ maximum  number  of  points  which  may  be  used  in 
each  initial  interval, 
mean  radius  in  log-normal  distribution 

Std.  deviation  in  log-normal  distribution  (as  input  data) ; 
used  later  as  ln(SIGMA) 

parameter  in  double  exponential  distrlb.,  and  in  exponential 
distributions 

parameter  in  several  distributions  and  in  exponential 
distributions 

parameter  in  several  distributions  and  in  exponential 
distributions 


:i 


j 


DENS 

DELRD 

VIS 

RC 

ALF 

GAM 

FOA 


RBARA 

RBARC 

DR(I) 

R(KK) 

F(KK) 

DELR 

DEN 

GNUM 

DENA 

DENC 

GNUMA 

GNUMC 

NRADII 

VOL 

VOLA 

VOLC 

AVOL 

OLSTAR 

OM2 

CTSUMT 

CSSUMT 

DENST 

CRSUMT 

EMM 

3H 

FH 

CH 

RE  XT 
KEXTT 
DRYVOL 
PSUHT  (J) 


particle  number  density 

increment  of  radius  for  Type  5 and  b distributions, 
visibility  in  km 

parameter  in  modified  gamma  distribution 
parameter  in  modified  gamma  distribution 
parameter  in  modified  gamma  distribution 

parameter  in  bimodal  log  normal  distribution  (number  density 
for  mode  "A") 

parameter  in  bimodal  log  normal  distribution  (number  density 
for  mode  "C")  ^ 

standard  deviation  user  supplied  bimodal  log  normal 
distribution 

standard  deviation  user  supplied  bimodal  log  normal 
distribution 

mean  radius  of  "accumulation' mode  for  bimodal  distributions 

mean  radius  of  "coarse"  mode  for  bimodal  distributions 

RR(I+1)  - RR(I)  ■ initial  interval  size,  1-1,  NI 

all  possible  values  of  radius,  KK-1,  NMAX 

value  of  distribution  function  at  each  R(KK) , KX“1,  NMAX 

increment  of  radius  for  arbitrary  distribution 

temporary  storage 

temporary  storage 

temporary  storage 

temporary  storage 

temporary  storage 

temporary  storage 

NMAX 

particle  volume 
particle  volume 
particle  volume 

analytically  evaluated  average  particle  volume 
value  of  expansion  coefficient  co^  for  total  aerosol 
value  of  expansion  coefficient  u for  total  aerosol 
value  of  extinction  cross-section  for  total  aerosol 
value  of  total  scattering  cross-section  for  total  aerosol 
value  of  particle  number  density  for  total  aerosol 
value  of  radar  cross-section  for  total  aerosol 
refractive  index  of  medium  in  which  scattering  particles 
are  deployed 

a factor  used  in  size  adjustment  at  non-zero  relative  humidity 
fractional  relative  humidity  (saturation  ratio) 
temporary  storage 

extinction  coefficient  for  a single  aerosol  component  (per  km) 
total  extinction  ceofficient  summed  over  all  components (per  km) 
total  dry  aerosol  particle  volume  in  (l*ns) 

final  average  phase  function  at  Jth  Gauss-Legendre  value  of  ;i 


The  pre-coded  (1DSTP  * 8,9,10)  values  of  SGA  and  SGC  are  the  natural 
logarithms  of  the  standard  deviations 


SYMBOL 


Explanation 


PC(J) 

RMS  (J) 

TEMK 

EMW 

CAYW 

RHOW 

EMA 

CAYA 

EMUA 

RHOA 

CONC 

DELTA 

VOL 

EM 

CAY 

FSUMG 

CTSUM 

CSSUM 

CRSUM 

OL1SUM 

OL2SUM 

PSUM(J) 

D 

RIT 

AC 

BC 

ALPHA 

QT 

QS 

QR 

P(J) 

PC(J) 

OISTAR 
0 2STAR 
FKK 
FKKA 

0 LGG.VOLHH 

oligg.olihh 
0L2GG,  OL2HH 
CTGG.CTHH 
CSGG.CSHH 


common  block  arrays  not  used  explicitly  In  subroutine  AGXPT2 
absolute  temperature  in  ®K 

real  part  of  refractive  index  of  water  at  temperature  TEMK 
imaginary  part  of  refractive  index  of  water  at  temperature  TEMK 
mass  density  (gm/cm3)  of  water  at  temperature  TEMK 
real  part  of  refractive  index  of  dry  aerosol  particle 
imaginary  part  of  refractive  index  of  dry  aerosol  particle 
(Hanel)  growth  factor  index  of  dry  aerosol  particle 
mass  density  (gm/cm3)  of  dry  aerosol  particle 
mass  concentration  (gm/cm3)  of  a dry  aerosol  component 
convergence  level  for  halving  method  of  integration 
total  volume  occupied  by  aerosol  material  distributed  in 
1 cm3  of  space 

real  part  of  refractive  index  of  actual  aerosol  particle 
imaginary  part  of  refractive  index  of  actual  aerosol  particle. 
Also,  CAY  - CAY/EM 

integral  of  particle  size  distribution  for  one  component 
of  an  aerosol 

extinction  cross-section  for  one  component  of  an  aerosol 
total  scattering  cross-section  for  one  component  of  an  aerosol 
radar  cross-section  for  one  component  of  an  aerosol 
value  of  £3  for  one  component  of  an  aerosol 
value  of  u>2  for  one  component  of  an  aerosol^ 
weighted  average  intensity  [ ( i^-t-ij)  / 2 ] at  J Gauss-Legendre 
value  of  UjJ  J»l,  IT. 
an  initial  interval  of  radius 

adjusted  radius,  according  to  relative  humidity 
parameter  in  formula  for  RIT 
parameter  in  formula  for  RIT 
size  parameter  ■ (2tt)  (radius) /(wavelength) 
extinction  efficiency  factor  for  a given  size  parameter 
scattering  efficiency  factor  for  a given  size  parameter 
radar  efficiency  factor  for  a given  size  parameter  ^ 
phase  function  factor  for  a given  size  parameter,  at  J value  of 
(phase  function) (distribution  function)  for  a given  size 
parameter,  at  Jth  value  of  y 
value  of  clj^  for  a given  size  parameter 
value  of  f°r  a Riven  size  parameter 
value  of  distribution  function  at  radius  R(KK) 

(FKK) (geometrical  cross-section  of  particle) 
partial  contributions  to  volume  of  average  particle,  for 
one  component  of  aerosol 

partial  contributions  to  average  average  particle,  for  one 
component  of  aerosol 

partial  contributions  to  average  average  particle,  for  one 
component  of  aerosol  \ 

partial  contributions  to  average  extinction  cross-section  for 
one  component  of  aerosol  \ 

partial  contributions  to  average  scattering  cross-section, 
for  one  component  of  aerosol 
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SYMBOL 


>lanation 


CRGG.CRHH 


FT 

VOLHHT 

OL1HHT 

OL2HHT 

CTHHT 

CRHHT 

CSHHT 


DENSC 

DENST 

PGG(J) 

PHH(J) 

PHHT(J) 

CONCT 

CATTN 

CATTNW 


partial  contributions  to  average  radar  cross-section,  for 
one  component  of  aerosol  gt 

fractional  change  in  contribution  to  volume  from  n to  n+1 
halving, for  one  initial  interval 
partial  contribution  to  integral  of  distribution  function 
contribution  to  volume  of  average  particle  from  one  original 
interval 

contribution  to  5^  of  average  particle  from  one  original  interval 
contribution  to  of  average  particle  from  one  original  interval 
contribution  to  average  extinction  cross-section  from  one 
original  interval 

contribution  to  average  radar  cross-section  from  one  original 
interval 

contribution  to  average  scattering  cross-section  from  one 
original  interval 

contribution  to  integral  of  distribution  function  from 
one  original  interval 

particle  number  density  for  a given  aerosol  component 
particle  number  density  for  total  aerosol 

partial  contributions  to  average  intensity  for  one 

aerosol  component  at  abscissa  value  y^ 

average  intensity  at  y^  for  one  original  interval 
total  aerosol  mass  concentration  (mg/m3),  summed  over  all 
components 

total  attenuation  coefficient  in  m2  per  mg  of  dry  aerosol 
material 

total  attenuation  coefficient  in  m2  per  mg  of  dry  aerosol  + 
accreted  water 


AGAUS10  (Main  Program)  - Simplified  Flow  Chart 


Read  Control  Parameter  Card 


Set-up  Size  Distribution 
Parameters 
(Subroutine  AGXPT1) 


Read  Wavelength,  Wavelength  Increment,  etc. 


Increment  Wavelength  Index  & ] EXIT 
TEST  FOR  EXIT 


Read  Optical  data,  etc. 
Perform  Mie  calculations  and 
integrations  over  sizes 

(Subroutine  AGXPT2) 


Calculate  Scattering  Fractions, 
etc.  at  given  wavelength  and  print 
results 

(Subroutine  AGXPT3) 


Sum  Results  over  Wavelength 


Return  for  next  wavelength 


Average  Results 
over 

Wavelength 


Print  Averaged 
Results 

(Subroutine  AGXPRT) 


Subroutine  AGXPT1  - Simplified  Flow  Chart 
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RETURN 


Subroutine  AGXPT2  - Simplified  Flowchart 
(Receives  control  data  and  size  distribution  data  from  main  program) 


Initialize  parameters  used  in  integrations 


Determine  optical  constants  and  mass  density 
for  liquid  water  at  given  wavelength 

(Subroutine  WATER;  skipped  for  IDSTP  • 12) 


Test  to  see  if  NINDX  aerosol  components 
have  been  treated  and  summed 


Read  optical  and  physical  parameter  card 
for  aerosol  (skipped  for  IDSTP  ■ 6) 


Do  final  renormalizf 
tion  over  components 


RETURN 


Begin  loop  over  size  intervals 


Treat  end-point  radii  for  size  interval  to 
get  first  estimate  for  Cext,  etc. 


Renormalize  results  over  component  NK 


(Loop  32) 
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A. 3.3  Input  Data  for  AGAUS9  and/or  AGAUS10 

The  data  and  control  cards  required  for  running  these  programs  are 
of  four  basic  types. 

TYPE  General  Description 

1 A set  of  Integers  which  select  certain  options  available 
within  the  program. 

2 A set  of  data  describing  the  parameters  of  the  particular 
size  distribution  function  to  be  used.  For  cases  other 
than  the  user-supplied  "arbitrary"  distribution  (IDSTP-O), 
only  one  card  of  this  type  is  needed. 

3 A set  of  data  describing  the  initial  wavelength  (vim)  to  be  used. 
Its  increment,  aerosol  number  density,  relative  humidity, 
atmospheric  temperature  (and,  for  AGAUS10,  IDSTP-6,  the 
desired  convergence  testing  level). 

4 A set  of  data  describing  the  optical  and  physical  properties 
of  the  aerosol  material.  If  the  aerosol  is  a mixture  of 
materials  of  unlike  properties,  more  than  one  card  of  this 
type  is  needed  for  each  card  of  type  3.  No  card  of  this 
type  is  used,  however,  for  runs  with  parameter  IDSTP 

(on  card  type  1)  is  equal  to  6. 

Remarks : 

1.  The  simplest  type  of  run  (IDSTP-6,  water  cloud  or  fog  model) 
requires  one  card  each  of  types  1,  2,  and  3.  For  other  IDSTP  choices, 
at  least  one  card  of  type  4 is  also  needed. 

2.  If  the  run  is  to  use  several  wavelengths,  then  at  least  one 
card  of  type  4 is  required  for  each  wavelength  (more  than  one  type  4 
card  will  be  needed  at  each  wavelength  if  the  aerosol  is  a multicomponent 
mixture) . 

3.  Two  special  "looping"  modes  of  AGADS9  and  AGAUS1Q  affect  the 
number  of  cards  of  type  3 which  are  needed: 

A)  For  runs  at  constant  relative  humidity  and  several 
wavelengths,  only  one  type  3 card  is  needed  or  permitted.  To  activate 
this  mode,  the  wavelength  increment  DWAVE  on  card  type  3 must  be  larger 
than  10~*. 

B)  Runs  at  constant  wavelength  and  a set  of  differing  values 
of  relative  humidity  require  one  set  of  cards  of  types  3 and  4 for  each 
value  of  relative  humidity.  This  option  is  invoked  by  setting  the 
parameter  "DWAVE"  equal  to  zero  on  the  first  and  all  subsequent  cards  of 
type  3. 
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Summary  of  Data  Card  Requirements 


CARD  TYPE  Input  Symbols 

1 IDSTT? , NRADI,  IDUMMY,  NWAVE,  NINDX,  MQRTE,  NCRDS,  NUNIT, 

IANG,  ISCAT  [FORMAT  (1115)] 

Explanations:  All  parameters  on  this  card  have  been 
defined  in  Section  A. 2.1,  with  the  exception  of 
NRADI.  NRADI  is  used  by  AGAUS10  only  for  the  case  IDSTP; 
its  valufe  for  other  IDSTP  choices  is  ignored 

2 Aerosol  size  distribution  data;  parameters  vary  with 
IDSTP  value;  see  Section  A. 2. 4. 

3 WAVE,  DWAVE,  RELHUM,  DENSH,  TEMP,  DELTAF 
[Format  (6D12.6)] 

Explanations : See  Sec.  A. 2.1  for  the  explanation  of  all 
quantities  other  than  DELTAF.  Parameter  DELTAF  is  the 
"convergence  level"  to  be  used  for  IDSTP-6  runs;  it's 
value  is  ignored  for  other  IDSTP  choices 

Cards  of  Type  3 may  be  repeated  (NWAVE- 1)  times  if  and 
only  if  DWAVE  < 1.0D-4 

4 EMA,  CAYA,  EMUA,  RHOA,  CONC,  DELTA 
[Format (4F10. 6,  2D15.7)] 

Explanations:  With  the  exception  of  DELTA,  all  parameters 
are  defined  in  Section  A. 2.5.  DELTA  is  the  desired 
(fractional)  convergence  level  (A);  it  should  be  positive. 

Remarks : 

A.  Card  Type  #4  is  not  required  and  is  never  read-in  if 
IDSTP^b  (water  cloud/fog  model)  since  the  relevant 
data  are  obtained  from  subroutine  WATER.  Since  card 
#4  is  not  used  for  that  case,  the  desired  value  of 
DELTA  must  be  inserted  on  card  type  #3  as  DELTAF. 

B.  In  the  case  IDSTP-12,  the  rain  model,  card  type  #4  must 
carry  the  optical  data  for  liquid  water  as  EMA  and  CAYA. 

The  reason  for  that  inconsistency  is  that  the  rain  model 
will  most  likely  be  usable  at  wavelengths  for  which  no 
data  appear  in  subroutine  WATER  (restricted  to  \ ^ 0.2  mm). 
[The  large  droplet  sizes  to  be  found  in  the  rain  model 
will  usually  cause  premature  failure  of  the  Mie  routine 

at  small  wavelengths]. 


4 
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C.  Card  Type  #4,  with  appropriately  adjusted  data 
values,  must  be  repeated  (NWAVE-1)  times. 


A. 3.4  Incidental  Remarks  Regarding  AGAUS10 

1)  In  AGAUS10,  the  average  "dry  volume  per  particle"  is  found  using 
all  available  values  for  particle  radii.  It  may  differ  from  that  printed 
in  routine  AGXPT2  since  that  routine  might  not  proceed  to  the  use  of  all 
available  values. 

2)  The  "volume"  convergence  tests  in  AGXPT2  use  the  volume  inferred 
after  any  growth  arising  from  non-zero  saturation  ratios  has  been  Included. 

3)  It  should  be  noted  that  the  convergence  tests  used  in  AGAUS10 
do  not  really  provide  tests  of  the  absolute  accuracies  achieved  for  the 
tested  quantity.  It  is  possible  for  "exit"  to  occur  (<$<A)  even  though  the 
use  of  "another"  halving  might  lead  to  <5>A  again.  If  runs  using  some  choice 
of  A should  lead  to  what  appear  to  be  "unusual"  or  unexpected  results, 

it  is  advised  khat,  as  a test,  the  case  be  re-run  using  a smaller  value 
of  A 

4)  Some  users  may  wish  to  explore  possible  increases  in  computation 
efficiency  which  might  result  from  changes  in  the  number  (NI)  of  size- 
intervals,  or,  the  ways  in  which  they  have  been  chosen  (in  AGXPT1). 


Comments  on  Usage  of  AGAUS9  and  AGAUS10  in  the  IANG-0  Mode 

AGAUS9  and  AGAUS10  have  been  coded  to  preserve  the  option  of  creating 
the  expansion  coefficients,  which  are  needed  by  the  ASL  multiple 
scattering  code  STAR04,  its  NHSU  version  AGSCAT,  and  the  NMSU/ASL  thermal 
emission  code  CLEM70.  Usage  of  the  coefficients  for  the  above  purposes, 
however,  requires  some  caution  on  the  part  of  the  user.  In  particular, 
users  must  not  assume  that  the  absence  of  abnormal  termination  of  runs  of 
AGAUS9  or  AGAUS10  guarantees  that  "all  is  well".  Past  experience  has  shown 
that  "warnings"  of  possible  problems  printed  by  those  programs  may  not 
be  noticed  or  taken  very  seriously.  In  an  attempt  to  overcome  those 
oversights,  additional  tests  and  warnings  have  been  incorporated  in  AGAUS9 
and  AGAUS10.  If  those  warnings  are  to  be  meaningful,  users  should 
examine  printed  outputs  carefully  for  such  warnings  before  using  the  results 
of  runs  as  Inputs  to  subsequent  codes. 


One  parcicular  point  which  should  be  checked  la  to  see  that  the 
quantity  printed  under  the  heading  ALBDO  agrees  reasonably  well  with  the 
zeroth  coefficient  (u>  ) of  the  expansion  coefficients.  Substantial 
disagreement  between  those  two  quantities  tsually  means  Chat  Che  paramecer 
"IT"  used  in  a run  was  "too  small"  to  acheive  a really  accurate 
reconstruction  of  the  phase  function  from  only  "it"  Legendre  expansion 
terms. 


A.  4 SAMPLES  OF  DATA  DECKS  FOR  AGAL'S9  AND  AGAUS10 
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Example  1:  IDSTP  - 0 (User  supplied  distribution  data) 


Column  1 


* 0 011 
0. 10000D-00 
0. 00000D-00 
0. 20000D-00 
0. 40000D-00 
0. 60000D-00 
0. 80000D-00 
1. 00000D-00 
0. 80000D-00 
0. 60000D-00 
0. 40000D-00 
0. 20000D-00 
0. 00000D-00 
10. 6000D+00 
1.9530 
10.6000D+00 
1.9520 
10. 6000D+00 
1.9530 
10. 6000D+00 
1.9530 
10. 6000D+00 
1.9530 
10. 6000D+00 
1.9530 
10. 6000D+00 
1.9530 


1 7 112345  0 0 05  01 

0. 02000D+00 


0. 10000D-05  00. OOOOD+OO  O.OOOOOD-OS  25. OOOOD+OO 
.468000  0.00000  1.87000  1.000000D-09  1. 0000000D-02 

0. 10000D-05  75. OOOOD+OO  0.00000D-05  25. OOOOD+OO 
.468000  0.15900  1.87000  1.000000D-09  1.0000000D-02 

0. 10000D-05  80. OOOOD+OO  O.OOOOOD-OS  25. OOOOD+OO 
.468000  0.15600  1.87000  1.000000D-09  1. 0000000D-02 

0.10000D-05  85. OOOOD+OO  0.00000D-05  25. OOOOD+OO 
.468000  0.15200  1.87000  1.000000D-09  1. 0000000D-02 

0. 10000D-05  90. OOOOD+OO  O.OOOOOD-OS  25.0000D-00 
.468000  0.14600  1.87000  1.000000D-09  1.0000000D-02 

0. 10000D-05  95. OOOOD+OO  O.OOOOOD-OS  25. OOOOD+OO 
.468000  0.14000  1.87000  1.000000D-09  1. OOOOOOOD-02 

0. 10000D-05  99. OOOOD+OO  O.OOOOOD-05  25. OOOOD+OO 
.468000  0.15800  1.87000  1.000000D-09  1. 0000000D-02 


This  example  is  set  up  to  read  11  distribution  data  cards, 
and  to  run  for  seven  values  of  relative  humidity  at  X - 10.6  uM. 
Since  MQRTE  - 12345,  detailed  Mie  results  for  each  of  the  11 
radii  will  be  printed.  Parameter  "it"  is  5,  and  IANG  - 1.  Note 
that  the  second  data  card  contains  RLO  (-0.1  uM)  and  DELLR  (-0.02). 
The  next  eleven  cards  contain  F(RLO) , F(RLO+DELLR,  etc.)  The 
remaining  cards  consist  of  seven  pairs  of  cards,  two  for  each  of 
the  seven  values  of  NWAVE.  In  all  cases,  DENSH  is  zero  which 
means  that  particle  number  densities  will  be  calculated  from 
RH0A  (1.87,  above),  CONC  (10~*  gm/cc,  above)  and  the  average 
particle  volume  computed  within  the  program. 
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Example  2:  IDSTP  ■ 1 (Log  normal  distribution) 


Column  S 
- 

1 100 
3. 700000-01 
10.  60000+00 
1.9530 
10. 60000+00 
1.9530 
10. 6000D+00 
1.9530 
10.60000+00 
1.9530 
10. 60000+00 
1.9530 
10.60000+00 
1.9530 
10.60000+00 
1.9530 


171  2345  0 0 05  00 

1.54000D+00  0. 00500D-00  1.00000D-00 
0. 10000D-05  OO.OOOOD+OO  O.OOOOOD-OS  25.0000D+00 
.468000  0.00000  1.87000  1.000000D-09  1.0000000D-02 

0.100000-05  75.00000+00  0.00000D-05  25.00000+00 
.468000  0.15900  1.87000  1.0000000-09  1.0000000D-02 

0. 10000D-05  80.00000+00  0.000000-05  25.0000D+00 
.468000  0.15600  1.87000  1.000000D-09  1.0000000D-02 

0. 100000-05  85.00000+00  0.000000-05  25.00000+00 
.468000  0.15200  1.87000  1.0000000-09  1. 0000000D-02 

0. 10000D-05  90.00000+00  0.00000D-05  25.00000+00 
.468000  0.14600  1.87000  1.000000D-09  1 . OOOOOOOD-02 

0.100000-05  95.00000+00  0.000000-05  25.00000+00 
.468000  0.14000  1.87000  1.0000000-09  1 . 0000000D-02 

0. 100000-05  99.00000+00  0.00000D-05  25.00000+00 
.468000  0.15800  1.87000  1.000000D-09  1.00000000-02 


Remarks:  This  run  also  uses  seven  values  of  relative  humidity. 

The  second  card  carries  RBAR,  SIGMA,  RL0  and  RHI.  If 
RLO  and  RHI  are  set  to  zero,  the  program  itself  will 
choose  RLO  and  RHI. 


Example  3:  IDSTP  - 5 (Modified  gamma  distribution) 

Column  5 

5 040  1 1 3 2345  0 0 01  00  0 

.0278520+00  1.225490+01  4.000000+00  .6000000+01  .1000000+01 
.7000000+00  . 000000D-04  .000000D+00  1.00000D+02  5.0000D+00 


1.  3300 

. 000000 

0.00000 

1.00000 

1.0000000+00 

1. OOOOOOOD-02 

1.4300 

. 000000 

0.00000 

1.00000 

1.000000D-01 

1. 0000000D-02 

1.5300 

. 050000 

0.00000 

1.00000 

1. 000OO0D-03 

1. 0000000D-02 

Remark:  This  example  sets  NINDX  • 3,  and  treats  a mixture  of 
three  aerosol  components  having  different  refractive 
indices  and  mass  concentrations. 
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Example  4:  IDSTP  ■ 6 (NMSU  water  cloud/fog  distribution) 


Column  5 
1 

6 000  1 5 1 2345  0 0 05  01  0 

. 027852D+00  1.22549D+01  4.00000D+00  .600000D+01  .100000D+01  1.00000D-06 
8 . OOOOOOD+OO  1. 00000D+00  .000000D+00  0. 00000D+O0  5.0000D+00  1.00000D-02 


Remark:  This  example  uses  the  same  distribution  parameters  as 
example  3,  but  also  carries  ELWC  - 10“*  gm/cm5  on  the 
second  data  card.  The  setup  is  for  runs  at  5 wavelengths, 
beginning  at  8.0  uM,  with  increments  of  1.0  yM. 


Example  5:  IDSTP  ■ 10  (Blmodal  distribution) 

Column  5 

1 

10  050 

10. 6000D+00 

1.9530 
10. 6000D+00 

1.9530 
10. 6000D+00 

1.9530 
10. 6000D+00 

1.9530 
10. 6000D+00 

1.9530 
10. 6000D+00 

1.9530 
10. 6000D+00 

1.9530 

Remarks:  This  run  set  up  is  similar  to  example  2,  except  that  it 
is  for  IDSTP  - 10. 


Example  6 IDSTP  ■ 12  (Marshall-Palmer  rain  model) 

Column  5 

r~ 

12  000  1 2 1 2345  0 0 05  01  0 

1. 00000D+O0 

1.000000D+03  1. 00000D+O3  .000000D+00  0. 00000D+O0  5.0000D+00  1.00000D-02 

1.3300  .000000  0.00000  1.00000  1.000000D-06  1. 0000000D-02 

1.3300  .000000  0.00000  1.00000  1.000000D-06  1. 0000000D-02 

Remarks:  The  second  card  carries  a rain  rate  of  1 mm  per  hour.  Two 
wavelengths  are  to  be  used  (1mm  and  2mm).  The  optical  data 
(last  two  cards)  are  merely  examples;  those  data  are  not 
really  correct  for  X - 1mm  and  2mm. 


. J 


171  2345  0 0 09  00 

. OOOOOOD-04  .000000D+00  .100000D-04  25.0000D+00 
.468000  0.00000  1.87000  1.000000D-09  1.0000000D-02 

0. 10000D-05  75. 0000D+00  0.00000D-05  25.0000D+00 
.468000  0.15900  1.87000  1.000000D-09  1. OOOOOOOD-02 

0. 10000D-05  80. 0000D+O0  0.00000D-05  25.0000D+00 
.468000  0.15600  1.87000  1.000000D-09  1. 0000000D-02 

0. 10000D-05  85.0000D+00  0.00000D-05  25.0000D+00 
.468000  0.15200  1.87000  1.000000D-09  1. 0000000D-02 

0. 10000D-05  90. 0000D+00  O.OOOOOD-05  25.0000D+00 
.468000  0.14600  1.87000  1.000000D-09  1. 0000000D-02 

0. 10000D-05  95. 0000D+O0  0.00000D-05  25.0000D+00 
.468000  0.14000  1.87000  1.000000D-09  1. 0000000D-02 

0. 10000D-05  99. 0000D+00  0.00000D-05  25.0000D+00 
.468000  0.15800  1.87000  1.000000D-09  1. 0000000D-02 
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B.l  PREPARATION  OF  TABLES 


Program  AGAUSLOhas  been  used  to  determine  the  scattering  fractions 
necessary  for  construction  of  tables.  Within  program  AGAUSX  the  following 
parameters  should  be  set  (additional  parameter  input  is  explained  in 
program  AGAUSX):  NCRDS  - 0,  IANG  - 1,  ISCAT  - 1,  IT  ■ 91.  The  output, 
comprised  of  the  scattering  fractions,  will  either  be  punched  on  cards 
or  written  on  NUNIT,  which  is  an  input  parameter  in  AGAUSX.  Once  the 
punched  deck  is  obtained  the  following  cards  must  precede  the  punched 
output: 

Card  1:  HEADER,  FORMAT  (20A4) , 

Card  2:  DWAVE,  FORMAT  (D12.6), 

Card  3:  CATTN,  FORMAT  (25x,D25. 14) , 

where  HEADER  is  for  informative  purposes,  DWAVE  is  the  wavelength  in 
microns  and  CATTN  is  the  attenuation  coefficient  in  m2  mg-1.  Additionally, 
a sentinel  card  must  be  inserted  at  the  end  of  the  deck  if  more  scattering 
fractions  (additional  wavelengths)  are  to  be  read  in.  The  sentinel 
card  is  (starting  in  column  1)  999  (no  decimal  point),  and  should  not 
follow  the  physically  last  data  card.  Program  CNTRL,  which  is  reproduced 
below,  will  then  create  a direct  access  (indexed)  file  which  may  be 
used  with  program  FIND.  A sample  run  stream  for  program  CNTRL  is 
provided  on  the  following  page.  Three  input  wavelengths  are  assumed 
in  the  example. 


B.2  DOCUMENTATION  FOR  PROGRAM  FIND 


Program  FIND  is  a stand-alone  program  or  may  easily  be  incorporated 
into  program  ASLSOM.  The  purpose  of  program  FIND  is  to  read  a given 
table  and  Interpolate  over  wavelength,  if  necessary,  to  find  precalculated 
values  of  the  scattering  fraction  and  attenuation  coefficient  (in  m2/mg). 
The  given  table  is  for  one  aerosol  model;  different  aerosols  require  a 
different  table.  Program  FIND  operates  by  using  a direct  access  (indexed) 
file  to  reduce  'look  up'  time.  The  control  procedure  is  as  follows: 
the  wavelengths  within  the  table  are  listed  in  increasing  order  and  are 
seventeen  records  long,  one  record  is  eighty  characters.  The  physical 
setup  of  one  'wavelength  block',  the  complete  information  pertaining  to 
the  wavelength  under  consideration,  contains  the  following  information: 

Record  1:  RLO,  RHI,  RC,  ALPHA,  GAMMA,  DENSITY 
Record  2:  WAVE1,  CATTN1,  FIRST 
Records  3-17:  SCATl(I) 
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Sample  Runstream  for  Program  CNTRL 


e FOR,  I 

•CNTRL 

CNTRL  DECK 

• 

• 

• 

@ MAP,  IN 

CNTL 

IN 

•CNTRL 

END 

<3  CAT,  P 

FILE. 

@ ASG,  A 

FILE. 

<3  USE  23.,  FILE. 

@ XQT 

CNTL 

DATA  DECK  1 

999 

DATA  DECK  2 


999 


DATA  DECK  3 


142 


01NMSU*APX(1)  . CNTRL 


POl  FOR  REWRITING  THE  OUTPUT  FROM  AGAUSX  INTO  ASLSOM  FORMAT 
FOR  TABLE  PREPARATION  TO  BE  USED  WITH  PROGRAM  FIND 
REAL *8  SCAT  (91) .CATTN, DWAVE 
DIMENSION  SSCAT  (91), HEADER  (20) 

N-0 

CHANGE  THE  FIRST  PARAMETER  IN  THE  DEFINE  FILE  23  STMT  TO  BE 
17  TIMES  THE  NUMBER  OF  INPUT  WAVELENGTHS  (in  this  case,  5 wavelengths) 
DEFINE  FILE  23(425, 80, L.IFIND) 
co:  TINUE 

READ  (5, 100, END-5)  (HEADER(I) ,1-1,20) 

WRITE  (23'IFIND,100)  (HEADER(I) ,1-1,20) 

READ  (5,102)  DWAVE 
WAVE-SNGL  (DWAVE) 

READ  (5,101)  CATTN 
CATN-SNGL(CATTN) 

READ  (5,102)  SCAT(l) 

SS CAT (l)-SNGL (SCAT (1)) 

WRITE  (23' IFIND.103)  WAVE,CATN,SSCAT(1) 

DO  2 1-2,91 

READ  (5,102)  SCAT(I) 

SSCAT(1)-SNGL(SCAT(I) ) 

DO  3 1-2,86,6 
K-I+5 

WRITE  (23'IFIND,103)  (SSCAT(J) , J-I ,K) 

READ  (5, 104, END-5)  NTEST 
N-N+l 

IF  (NTEST. EQ. 999)  GO  TO  6 
PRINT  104, N 
3 FORMAT  (20A4) 

1 FORMAT  (25X.D25.14) 

2 FORMAT  (D12.6) 

3 FORMAT  (6(E12.6,1X)) 

4 FORMAT  (IX,  'END  PGM  CHANGE' ,2X, 13, 

+ 'FILES  HAVE  BEEN  concatenated') 
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where 


RLO  ■ minimum  aerosol  radius  (pm) 

RHI  “ maximum  aerosol  radius  (jjm) 

RC  ■ mode  aerosol  radius  (van) 

ALPHA  ■ parameter  in  the  modified  gamma  distribution 
GAMMA  ■ parameter  in  the  modified  gamma  distribution 
DENSITY  - aerosol  density  (If cm”3) 

WAVE1  ■ wavelength  under  consideration  (ym) 

CATTN1  - attenuation  coefficient  (m2  mg-1) 

FIRST  - scattering  fraction  at  zero  degrees 
SCATl(I)  - scattering  fraction  at  2°  to  180°  in  2°  increments. 

The  formats  are  as  follows: 

* 

Record  1 : 6E12.6 
Record  2 : 3(E12.6,lx) 

Records  3-17  : 6(E12.6,lx) 


record  1 is  read  as  Header  [Format  (A4) ] and  is  not  currently  used. 

Thus  given  the  proper  input,  to  be  subsequently  explained,  program 
FIND  will  either;  1)  find  the  input  wavelength  to  ± 2%  or  2)  find  the 
next  largest  wavelength  than  the  input  wavelength  and  then  interpolate 
over  the  internal  wavelength  range.  In  either  case  the  output  is 
either  punched  or  written  as  an  auxiliary  unit  under  the  following 
format  (ASLSOM  format): 

Records  1-91:  scattering  fractions,  FORMAT  (E12.6) 

Record  92  : attenuation  coefficient,  FORMAT  (10x,E10.4) 

Room  has  been  left  on  the  92n<1  record  to  insert  additional  (needed) 
parameters  in  their  proper  format  without  having  to  relocate  the  physical 
position  of  the  attenuation  coefficient. 

The  input  for  program  FIND  is  as  follows: 

one  card:  WAVE,  NUNIT:  FORMAT  (F7.3.I2) 

where  WAVE  - wavelength  (pm)  for  which  scattering  fractions  and 

attenuation  coefficient  are  desired 
NUNIT  ■ auxiliary  unit  for  scattering  fractions  and  attenuation 
coefficient  to  be  written  upon;  the  default  value  of 
NUNIT  (zero  or  blank)  is  the  card  punch  (unit  4 on 
Univac  1100  series). 

In  order  to  run  program  FIND  two  additional  commands  are  needed. 

One  must  assign  one  of  the  two  existing  fog  model  files,  8201NMSU*WSMRF01 
or  8201NMSU*WSMRF02  (see  previous  explanation  for  the  difference  in  fog 
models)  and  then  attach  a use  assigned  name  of  23  to  the  particular  file 
under  consideration.  Therefore  a typical  run  stream  might  appear  as 
follows: 


@ASG,  A 8201NMSU*WSMRF01 . 

0USE  23.,  8201NMSU*WSMRF01. 

<§FOR,  I .FIND 
CARD  DECK 

@MAP,  IN  FND 
IN  .FIND 

END 

@XQT  FND 

+01.06022 

@FIN 

This  would  assign  fog  model  WSMRF01  and  write  the  scattering  fractions 
and  attenuation  coefficient  for  1.06  urn  on  unit  22.  If  the  input  card 
had  read  +01.060  the  output  would  have  been  punched  on  cards. 


j 

I 
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Appendix  C 


FORMAL  SOLUTION  OF  THE  RADIATIVE  TRANSFER  EQUATION 


( y\^Z 

/ /.s'  , • —y 

* / 


It  is  desired  to  show  that  Eqs.  (3), 
(4)  are  solutions  of  Eqs.  (2),  (1), 
respectively.  In  order  to  do  this,  it 
is  first  necessary  to  show  that 


(C-l) 


where  £(r,s)  is  the  distance  from  the 
point  r In~the  cloud  back  along  (-s) 
to  the^poinu  R on  the  bounding  surface  E, 
as  illustrated  in  Fig.  Bl.  Note  that 


R ■ r - s i(r ,s) . 


(C-2) 


Let  the  surface  V be  given  by  f(R)  ■ 0 • f(r-is),  f(R+dR)  - 0 - 


f (R)  + dRi(9f/3R1) , so 


0 - drk(3Ri/3rk)(3f/3R1)  - drk(5kl-3ki  s^Of/SR^, 


(C-3) 


where  summation  rotation  is  used  (repeated  Indices  summed  over).  Since 
the  drk  are  arbitrary,  it  must  be  true  that  (5  -(3  £.)s  )(3f/3R  ) ■ 0. 

But  the  3f/3R^  are  not  zero;  therefore  the  determinant  of  the  matrix 
A^  ■ 6kj-3kS.  si  must  vanish.  But 


det  A ■ 1 - s.3,£  ■ 1 - s*V  - 0, 
i i ~ ~r  ’ 


(C-4) 


which  completes  the  desired  proof.  As  an  example  of  the  relation 
s*V  4 - 1,  consider  a plane  parallel  cloud,  as  in  Fig.  B2.  Then 

, i 4(r,s)  - z/cos6  ■ z/u 


s ■ i sin0  cos<t>  + J sin6  sin<J>  + ky 


V^Jl  ■ k 34/3z  ■ ky  . 


9*Vr<,  - kyky-1  - 1. 

Here,  (i,j,k)  are  the  usual  orthogonal 
triplet^of  cartesian  basis  vectors. 
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Now,  consider  the  expression  for  the  incident  beam  intensity  in  the 

•v  it  K 


®*!rrin  " ~Y(r!r4)Iin  + <exP-Y*>  (3lin <5 » £> / 3Ric)  (C-5) 

bnt  s-7.«,  - 1,  and  st  3^/3^  - 3^)  - sk(l-s-7rl)  - 0, 

so  indeed  Iln(r,s)  satisfies  Eq.  (2). 

Now,  consider  Eq.  (4),  the  expression  for  the  total  intensity.  Since 
s is  a unit  vector,  it  is  clear  that 


s*7rI(r-Cs,s')  - -3l(r-^s,s’)/3C.  (C-6) 

Using  this,  and  Eq.  (4),  and  s*7r«,(r,s)  - 1, 

s*7rr(r,s)  - s-7rIln  + (4ir)-‘y  e'Y*  |dfl#l  p(s,s')I(r-Jls,s') 

l 

- (4ir)-lYJd^  |^[e"Y5  jdflg,  p(ss')I(r-£s,s')  ] 

o 

- (Air)-1  y2Jd£  e~YC  jd^.pCs.s'JKr-Cs.s') . 

o 

The  last  term  on  the  RHS  of  this  equation  is  just  -y(I(r,s)-I  (r  s))- 
the  third  term  on  the  RHS  is  integrable,  and  yields  two'terms^oSe'of  ’ 
which  cancels  the  second  term.  These  results  yield 

s-^Kr.s)  - ■•Vn+(4ir)"lYjdn.'P<;»;,>I^'!,)-Y(I(r.;)-Iln(rt:)). 
But  s‘7^1^  + ylln  ■ 0,  (Eq.  2),  so  the  final  result  is 

r!rI(E*!)  +YKr,s)  - (4ir)-‘yfdn  .pCs.s'JKr.s'),  (C-7) 

J 9 — m, 

which  is  just  Eq.  (1). 
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The  formal  implicit  solution  of  Eq.  (C-l)  given  by  Eq.  (C-4)  is 
simply  a coordinate-dependent  realization  of  a familiar  form  given  by 
Chandrasekhar  [32].  The  authors  of  this  report  believe  that  this  form 
(Eq.  C-4)  has  not  previously  been  reported  in  the  literature. 


Addendum:  Most  routines  used  and  developed  under  this  contract  have 
been  placed  in  a special  file  in  WSMR's  system  B.  That 
file  is  "8201NMSU*NMSFINAL" . Some  notes  on  the  contents 
of  that  file  follow. 


Notes  on  File  8201NMSU*NMSFINAL 

This  file  contains  various  source  codes,  data  elements,  mapping 
and  running  elements  developed  or  modified  during  1978. 

There  are  two  principal  types  of  programs: 

(1)  ATRAN  - the  high  resolution  code 

(2)  assorted  single-scattering  (Mie)  codes. 

ATRAN:  The  main  program  is  ATRAN14.  A sample  1108  run-stream  will 
be  found  in  element  ATF.UNF.  That  element  either  contains  (or  will  direct 
one  to  an  element  which  contains)  compilation,  mapping  and  running 
instructions.  Reference  to  ATRUNF  or  elements  called  by  it  (via  @ADD 
statements)  will  define  which  source  codes,  etc.,  in  file  8201NMSU*NMSFINAL 
belong  to  "ATRAN". 

[Requires  about  63K-words  of  core! 

Mie-Codes : There  are  several  of  those: 

1)  AGAUS9  - an  NMSU  version  of  PGAUSS,  including  computation 

of  "scattering  fractions"  (documented  in  final  report). 

Uses  subroutines : AG9PT1 , AG9PT2 , AG9PT3 , AG9PRT , GAUS , GUSF.T , 

WATER, MIEGX, TIMB .ANGLE . 

Runstream:  element  * RUNAG9 

Sample  Data:  elements  - TYPE0,TYPE1,TYPE6,TYPE10,TYPE12,AGXC1 
(type  - 5) . 

[Requires  about  60K  of  core] 

2)  AGAUSX  (documented  in  the  final  report  as  AGAUS10)  - a "halving" 

version  of  AGAUS9. 

Uses  Subroutines:  AGXPT1 , AGXPT2 , AGXPT3 , AGXPRT , GAUS , GUSET .ANGLE , 
MIEGX, TIMB .WATER 
Sample  Runstream:  element  ■ RUNAGX 

[Requires  about  60K  of  core] 

Sample  Data  Elements:  Same  as  for  AGAUS9 


3)  MAIN1S  (AGAUS10S,  no  documentation  furnished).  This 

Is  a version  of  AGAUS10  which  uses  double-precision  arithmetic 
only  in  the  Mie-subroutine  (MIEGXS)  and  the  routine  (SGUSET) 
used  to  get  the  Gauss-Legendre  quadrature  weights  and 
abscissae.  It  is  logically  identical  to  AGAUSX,  except 
that  it  performs  the  convergence  checks  on  the  ave.  Intensity 
as  well  as  on  the  extinction  coefficient.  For  similar 
conditions  AGAUS10S  results  will  agree  with  those  of  AGAUS10 
to  six  or  more  digits,  but  it  needs  only  about  40K  words  of  core. 

Uses  Subroutines : PARTIS , PART2S , PART3S , AGP RTS .MIEGXS , SGAUS , 

SGUSET .ANGLES , TIMBS , WATERS 
Sample  Rustream;  element  ■ RUN10S 
Sample  Data:  same  as  AGAUS9  or  AGAUS10 

4)  AGAUSXL  - A version  of  AGAUS10  which  uses  the  RRA  Mie-routine 
instead  of  Querfeld's  forward-recursion  variant  (MIEGX) . The 
RRA  routine  DOWN42  is  used  in  a slightly  modified  form  called 
RRAMIE. 

Requires  Subroutines:  AGXPT1,AGXPT2L,AGXPT3,AGXPRT,GAUS , 

GUSET , WATER , TIMB , RRAMIE .SETUP , ANGLE 
Sample  Runstream:  RUNXL 
Sample  Data:  Same  as  AGAUS9.10.10S. 

[Core  Requirement  s 64K] 

Note:  No  documentation  has  been  prepared  for  AGAUSXL,  but 

most  of  the  logic  is  the  same  as  for  AGAUS10.  However, 
some  improvements  contained  in  AGAUS10  may  not  have 
been  incorporated  into  AGAUSXL  and  subroutine  AGXPT2L. 

5)  MIE2  - this  is  the  RRA  aerosol  extinction  routine  supplied  to 
NMSU  by  ASL.  It  is  in  file  8201NMSU*NMSFINAL  more  for  reference 
than  for  any  other  reason.  It  uses  elements  FIRST(=  DOWN42) 
and  SECOND (=  IPLAW) . A sample  data  element  is  ME2C1. 


^SPECIAL  NOTES# 

Programs  ATRAN14  and  AGAUS10S  exist  at  NMSU  in  IBM  360/370  compatible 
versions.  The  other  routines  listed  above  are  extant  only  in  Univac  1108 
tested  forms. 

[The  Univac  1108  version  of  ATRAN  at  WSMR  expects  a 7-track  (BCD) 
copy  of  the  AFGL  data— tape,  while  the  IBM  360/370  version  uses  a 9-track 
(EBCDIC)  version  of  the  tape.  Some  Job  Control  Language  (JCL)  differences 
also  exist,  but  IBM  360/370  JCL  instructions  could  be  supplied  if  desired]. 

AGAUS10S  has  also  been  successfully  modified  at  NMSU  to  run  (with 
slightly  restricted  capabilities)  on  a DEC/PDP-1160  computer.  That  version 
will  run  within  a 32K-word  memory  allocation.  Program  listings  could  be 
supplied  to  ASL  without  much  difficulty,  but  tape  or  card  copies  are  not 
presently  available. 


